Changeset 7522


Ignore:
Timestamp:
2017-01-02T11:06:49+01:00 (4 years ago)
Author:
cetlod
Message:

3.6 stable : update the offline routines to be able to run passive tracers offline with linear free surface, see ticket #1775

Location:
branches/2015/nemo_v3_6_STABLE/NEMOGCM
Files:
11 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/ARCH/arch-macport_osx.fcm

    r6204 r7522  
    5454%CPP               cpp-mp-4.8 
    5555%FC                mpif90  
    56 %FCFLAGS             -fdefault-real-8 -O3 -funroll-all-loops -fcray-pointer  
     56%FCFLAGS             -fdefault-real-8 -O3 -funroll-all-loops -fcray-pointer -ffree-line-length-none 
    5757%FFLAGS              %FCFLAGS 
    5858%LD                  %FC 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist_cfg

    r5385 r7522  
    101101!          !  file name  ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask ! 
    102102!          !             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  ! filename      ! 
    103    sn_tem  = 'dyna_grid_T' ,    120            , 'votemper' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
    104    sn_sal  = 'dyna_grid_T' ,    120            , 'vosaline' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
    105    sn_mld  = 'dyna_grid_T' ,    120            , 'somixhgt' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
    106    sn_emp  = 'dyna_grid_T' ,    120            , 'sowaflup' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
    107    sn_fmf  = 'dyna_grid_T' ,    120            , 'iowaflup' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
    108    sn_ice  = 'dyna_grid_T' ,    120            , 'soicecov' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
    109    sn_qsr  = 'dyna_grid_T' ,    120            , 'soshfldo' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
    110    sn_wnd  = 'dyna_grid_T' ,    120            , 'sowindsp' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
    111    sn_uwd  = 'dyna_grid_U' ,    120            , 'vozocrtx' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
    112    sn_vwd  = 'dyna_grid_V' ,    120            , 'vomecrty' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
    113    sn_wwd  = 'dyna_grid_W' ,    120            , 'vovecrtz' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
    114    sn_avt  = 'dyna_grid_W' ,    120            , 'voddmavs' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
    115    sn_ubl  = 'dyna_grid_U' ,    120            , 'sobblcox' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
    116    sn_vbl  = 'dyna_grid_V' ,    120            , 'sobblcoy' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
    117    sn_ahu  = 'dyna_grid_U' ,    120            , 'vozoahtu' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
    118    sn_ahv  = 'dyna_grid_V' ,    120            , 'vomeahtv' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
    119    sn_ahw  = 'dyna_grid_W' ,    120            , 'voveahtz' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
    120    sn_eiu  = 'dyna_grid_U' ,    120            , 'vozoaeiu' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
    121    sn_eiv  = 'dyna_grid_V' ,    120            , 'vomeaeiv' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
    122    sn_eiw  = 'dyna_grid_W' ,    120            , 'soleaeiw' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
     103   sn_tem  = 'dyna_grid_T' ,    120            , 'votemper'   ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
     104   sn_sal  = 'dyna_grid_T' ,    120            , 'vosaline'   ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
     105   sn_mld  = 'dyna_grid_T' ,    120            , 'somixhgt'   ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
     106   sn_emp  = 'dyna_grid_T' ,    120            , 'sowaflup'   ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
     107   sn_fmf  = 'dyna_grid_T' ,    120            , 'iowaflup'   ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
     108   sn_ice  = 'dyna_grid_T' ,    120            , 'soicecov'   ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
     109   sn_qsr  = 'dyna_grid_T' ,    120            , 'soshfldo'   ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
     110   sn_wnd  = 'dyna_grid_T' ,    120            , 'sowindsp'   ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
     111   sn_uwd  = 'dyna_grid_U' ,    120            , 'uocetr_eff' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
     112   sn_vwd  = 'dyna_grid_V' ,    120            , 'vocetr_eff' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
     113   sn_wwd  = 'dyna_grid_W' ,    120            , 'wocetr_eff' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
     114   sn_avt  = 'dyna_grid_W' ,    120            , 'voddmavs'   ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
     115   sn_ubl  = 'dyna_grid_U' ,    120            , 'sobblcox'   ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
     116   sn_vbl  = 'dyna_grid_V' ,    120            , 'sobblcoy'   ,  .true.    , .true. ,   'yearly'  , ''       , ''    , '' 
    123117! 
    124    cn_dir      = './'      !  root directory for the location of the dynamical files 
    125    ln_degrad   =  .false.  !  flag for degradation -                requires ("key_degrad") 
    126    ln_dynwzv   =  .true.   !  computation of vertical velocity instead of using the one read in file 
    127    ln_dynbbl   =  .true.   !  bbl coef are in files, so read them - requires ("key_trabbl") 
     118   cn_dir          = './'       !  root directory for the location of the dynamical files 
     119   ln_dynrnf       =  .false.   !  runoffs option enabled (T) or not (F) 
     120   ln_dynrnf_depth =  .false.   ! runoffs is spread in vertical (T) or not (F) 
     121!   fwbcorr      = 3.786e-06    ! annual global mean of empmr for ssh correction 
    128122/ 
    129123!----------------------------------------------------------------------- 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OFF_SRC/domrea.F90

    r5504 r7522  
    1717   USE lib_mpp         ! distributed memory computing library 
    1818 
     19   USE iom 
    1920   USE domstp          ! domain: set the time-step 
    2021 
     
    7374 
    7475      CALL dom_nam      ! read namelist ( namrun, namdom, namcla ) 
     76      CALL dom_msk      ! Masks 
     77      CALL dom_hgr      ! Horizontal grid 
    7578      CALL dom_zgr      ! Vertical mesh and bathymetry option 
    76       CALL dom_grd      ! Create a domain file 
    77  
    78      ! 
    79       ! - ML - Used in dom_vvl_sf_nxt and lateral diffusion routines 
    80       !        but could be usefull in many other routines 
     79      ! 
    8180      e12t    (:,:) = e1t(:,:) * e2t(:,:) 
    8281      e1e2t   (:,:) = e1t(:,:) * e2t(:,:) 
     
    9190      re1v_e2v(:,:) = e1v(:,:) / e2v(:,:) 
    9291      ! 
    93       hu(:,:) = 0._wp                          ! Ocean depth at U- and V-points 
    94       hv(:,:) = 0._wp 
    95       DO jk = 1, jpk 
    96          hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk) 
    97          hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk) 
    98       END DO 
    99       !                                        ! Inverse of the local depth 
    100       hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask(:,:,1) ) * umask(:,:,1) 
    101       hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask(:,:,1) ) * vmask(:,:,1) 
    102  
    10392      CALL dom_stp      ! Time step 
    104       CALL dom_msk      ! Masks 
    10593      CALL dom_ctl      ! Domain control 
    10694 
     
    178166      nstocklist = nn_stocklist 
    179167      nwrite = nn_write 
    180  
    181  
    182168      !                             ! control of output frequency 
    183169      IF ( nstock == 0 .OR. nstock > nitend ) THEN 
     
    222208904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp ) 
    223209      IF(lwm) WRITE ( numond, namdom ) 
     210 
    224211 
    225212      IF(lwp) THEN 
     
    321308   END SUBROUTINE dom_nam 
    322309 
     310   SUBROUTINE dom_msk 
     311      !!--------------------------------------------------------------------- 
     312      !!                 ***  ROUTINE dom_msk  *** 
     313      !! ** Purpose :  Read the NetCDF file(s) which contain(s) all the 
     314      !!      ocean mask informations and defines the interior domain T-mask. 
     315      !! 
     316      !! ** Method  :  Read in a file all the arrays generated in routines 
     317      !!               dommsk:   'mask.nc' file 
     318      !!              The interior ocean/land mask is computed from tmask 
     319      !!              setting to zero the duplicated row and lines due to 
     320      !!              MPP exchange halos, est-west cyclic and north fold 
     321      !!              boundary conditions. 
     322      !! 
     323      !! ** Action :   tmask_i  : interiorland/ocean mask at t-point 
     324      !!               tpol     : ??? 
     325      !!---------------------------------------------------------------------- 
     326      ! 
     327      INTEGER  ::  inum   ! local integers 
     328      INTEGER  ::   ji, jj, jk                   ! dummy loop indices 
     329      INTEGER  ::   iif, iil, ijf, ijl       ! local integers 
     330      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk 
     331      ! 
     332      !!--------------------------------------------------------------------- 
     333       
     334 
     335 
     336      IF(lwp) WRITE(numout,*) 
     337      IF(lwp) WRITE(numout,*) 'dom_rea : read NetCDF mesh and mask information file(s)' 
     338      IF(lwp) WRITE(numout,*) '~~~~~~~' 
     339 
     340      CALL wrk_alloc( jpi, jpj, zmbk ) 
     341      zmbk(:,:) = 0._wp 
     342 
     343      IF(lwp) WRITE(numout,*) '          one file in "mesh_mask.nc" ' 
     344      CALL iom_open( 'mask', inum ) 
     345 
     346         !                                                         ! masks (inum2)  
     347      CALL iom_get( inum, jpdom_data, 'tmask', tmask ) 
     348      CALL iom_get( inum, jpdom_data, 'umask', umask ) 
     349      CALL iom_get( inum, jpdom_data, 'vmask', vmask ) 
     350      CALL iom_get( inum, jpdom_data, 'fmask', fmask ) 
     351 
     352      CALL lbc_lnk( tmask, 'T', 1._wp )    ! Lateral boundary conditions 
     353      CALL lbc_lnk( umask, 'U', 1._wp )       
     354      CALL lbc_lnk( vmask, 'V', 1._wp ) 
     355      CALL lbc_lnk( fmask, 'F', 1._wp ) 
     356 
     357#if defined key_c1d 
     358      ! set umask and vmask equal tmask in 1D configuration 
     359      IF(lwp) WRITE(numout,*) 
     360      IF(lwp) WRITE(numout,*) '**********  1D configuration : set umask and vmask equal tmask ********' 
     361      IF(lwp) WRITE(numout,*) '**********                                                     ********' 
     362 
     363      umask(:,:,:) = tmask(:,:,:) 
     364      vmask(:,:,:) = tmask(:,:,:) 
     365#endif 
     366 
     367#if defined key_degrad 
     368      CALL iom_get( inum, jpdom_data, 'facvolt', facvol ) 
     369#endif 
     370 
     371      CALL iom_get( inum, jpdom_data, 'mbathy', zmbk )              ! number of ocean t-points 
     372      mbathy (:,:) = INT( zmbk(:,:) ) 
     373      misfdep(:,:) = 1                                               ! ice shelf case not yet done 
     374       
     375      CALL zgr_bot_level                                             ! mbk. arrays (deepest ocean t-, u- & v-points 
     376 
     377      !                                     ! ============================ 
     378      !                                     !        close the files  
     379      !                                     ! ============================ 
     380 
     381      ! 
     382      ! Interior domain mask (used for global sum) 
     383      ! -------------------- 
     384      ssmask(:,:)  = tmask(:,:,1) 
     385      tmask_i(:,:) = tmask(:,:,1) 
     386      iif = jpreci                        ! thickness of exchange halos in i-axis 
     387      iil = nlci - jpreci + 1 
     388      ijf = jprecj                        ! thickness of exchange halos in j-axis 
     389      ijl = nlcj - jprecj + 1 
     390      ! 
     391      tmask_i( 1 :iif,   :   ) = 0._wp    ! first columns 
     392      tmask_i(iil:jpi,   :   ) = 0._wp    ! last  columns (including mpp extra columns) 
     393      tmask_i(   :   , 1 :ijf) = 0._wp    ! first rows 
     394      tmask_i(   :   ,ijl:jpj) = 0._wp    ! last  rows (including mpp extra rows) 
     395      ! 
     396      !                                   ! north fold mask 
     397      tpol(1:jpiglo) = 1._wp 
     398      !                                 
     399      IF( jperio == 3 .OR. jperio == 4 )   tpol(jpiglo/2+1:jpiglo) = 0._wp    ! T-point pivot 
     400      IF( jperio == 5 .OR. jperio == 6 )   tpol(     1    :jpiglo) = 0._wp    ! F-point pivot 
     401      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot: only half of the nlcj-1 row 
     402         IF( mjg(ijl-1) == jpjglo-1 ) THEN 
     403            DO ji = iif+1, iil-1 
     404               tmask_i(ji,ijl-1) = tmask_i(ji,ijl-1) * tpol(mig(ji)) 
     405            END DO 
     406         ENDIF 
     407      ENDIF  
     408      ! 
     409      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at 
     410      ! least 1 wet u point 
     411      DO jj = 1, jpjm1 
     412         DO ji = 1, fs_jpim1   ! vector loop 
     413            umask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:))) 
     414            vmask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:))) 
     415         END DO 
     416         DO ji = 1, jpim1      ! NO vector opt. 
     417            fmask_i(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   & 
     418               &            * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:))) 
     419         END DO 
     420      END DO 
     421      CALL lbc_lnk( umask_i, 'U', 1._wp )      ! Lateral boundary conditions 
     422      CALL lbc_lnk( vmask_i, 'V', 1._wp ) 
     423      CALL lbc_lnk( fmask_i, 'F', 1._wp ) 
     424 
     425      ! 3. Ocean/land mask at wu-, wv- and w points  
     426      !---------------------------------------------- 
     427      wmask (:,:,1) = tmask(:,:,1) ! ???????? 
     428      wumask(:,:,1) = umask(:,:,1) ! ???????? 
     429      wvmask(:,:,1) = vmask(:,:,1) ! ???????? 
     430      DO jk = 2, jpk 
     431         wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1) 
     432         wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1)    
     433         wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1) 
     434      END DO 
     435      ! 
     436      CALL wrk_dealloc( jpi, jpj, zmbk ) 
     437      ! 
     438      CALL iom_close( inum ) 
     439      ! 
     440   END SUBROUTINE dom_msk 
     441 
     442   SUBROUTINE zgr_bot_level 
     443      !!---------------------------------------------------------------------- 
     444      !!                    ***  ROUTINE zgr_bot_level  *** 
     445      !! 
     446      !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays) 
     447      !! 
     448      !! ** Method  :   computes from mbathy with a minimum value of 1 over land 
     449      !! 
     450      !! ** Action  :   mbkt, mbku, mbkv :   vertical indices of the deeptest 
     451      !!                                     ocean level at t-, u- & v-points 
     452      !!                                     (min value = 1 over land) 
     453      !!---------------------------------------------------------------------- 
     454      ! 
     455      INTEGER ::   ji, jj   ! dummy loop indices 
     456      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk 
     457      !!---------------------------------------------------------------------- 
     458 
     459      ! 
     460      IF(lwp) WRITE(numout,*) 
     461      IF(lwp) WRITE(numout,*) '    zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels ' 
     462      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~' 
     463      ! 
     464      CALL wrk_alloc( jpi, jpj, zmbk ) 
     465      ! 
     466      mbkt(:,:) = MAX( mbathy(:,:) , 1 )    ! bottom k-index of T-level (=1 over land) 
     467      mikt(:,:) = 1 ; miku(:,:) = 1; mikv(:,:) = 1; ! top k-index of T-level (=1 over open ocean; >1 beneath ice shelf) 
     468      !                                     ! bottom k-index of W-level = mbkt+1 
     469      DO jj = 1, jpjm1                      ! bottom k-index of u- (v-) level 
     470         DO ji = 1, jpim1 
     471            mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  ) 
     472            mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  ) 
     473         END DO 
     474      END DO 
     475      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk  
     476      zmbk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 
     477      zmbk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 
     478      ! 
     479      CALL wrk_dealloc( jpi, jpj, zmbk ) 
     480      ! 
     481   END SUBROUTINE zgr_bot_level 
     482 
     483   SUBROUTINE dom_hgr 
     484      !!---------------------------------------------------------------------- 
     485      !!                  ***  ROUTINE dom_hgr  *** 
     486      !!                    
     487      !! ** Purpose :  Read the NetCDF file(s) which contain(s) all the 
     488      !!      ocean horizontal mesh informations  
     489      !! 
     490      !! ** Method  :   Read in a file all the arrays generated in routines 
     491      !!                domhgr:   'mesh_hgr.nc' file 
     492      !!---------------------------------------------------------------------- 
     493      !! 
     494      INTEGER ::   ji, jj   ! dummy loop indices 
     495      INTEGER  ::  inum    ! local integers 
     496      !!---------------------------------------------------------------------- 
     497 
     498      IF(lwp) WRITE(numout,*) 
     499      IF(lwp) WRITE(numout,*) 'dom_grd_hgr : read NetCDF mesh and mask information file(s)' 
     500      IF(lwp) WRITE(numout,*) '~~~~~~~' 
     501 
     502      IF(lwp) WRITE(numout,*) '          one file in "mesh_mask.nc" ' 
     503      CALL iom_open( 'mesh_hgr', inum ) 
     504 
     505      !                                                         ! horizontal mesh (inum3) 
     506      CALL iom_get( inum, jpdom_data, 'glamt', glamt ) 
     507      CALL iom_get( inum, jpdom_data, 'glamu', glamu ) 
     508      CALL iom_get( inum, jpdom_data, 'glamv', glamv ) 
     509      CALL iom_get( inum, jpdom_data, 'glamf', glamf ) 
     510 
     511      CALL iom_get( inum, jpdom_data, 'gphit', gphit ) 
     512      CALL iom_get( inum, jpdom_data, 'gphiu', gphiu ) 
     513      CALL iom_get( inum, jpdom_data, 'gphiv', gphiv ) 
     514      CALL iom_get( inum, jpdom_data, 'gphif', gphif ) 
     515 
     516      CALL iom_get( inum, jpdom_data, 'e1t', e1t ) 
     517      CALL iom_get( inum, jpdom_data, 'e1u', e1u ) 
     518      CALL iom_get( inum, jpdom_data, 'e1v', e1v ) 
     519       
     520      CALL iom_get( inum, jpdom_data, 'e2t', e2t ) 
     521      CALL iom_get( inum, jpdom_data, 'e2u', e2u ) 
     522      CALL iom_get( inum, jpdom_data, 'e2v', e2v ) 
     523 
     524      CALL iom_get( inum, jpdom_data, 'ff', ff ) 
     525 
     526 
     527      ! Control printing : Grid informations (if not restart) 
     528      ! ---------------- 
     529 
     530      IF(lwp .AND. .NOT.ln_rstart ) THEN 
     531         WRITE(numout,*) 
     532         WRITE(numout,*) '          longitude and e1 scale factors' 
     533         WRITE(numout,*) '          ------------------------------' 
     534         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   & 
     535            glamv(ji,1), glamf(ji,1),   & 
     536            e1t(ji,1), e1u(ji,1),   & 
     537            e1v(ji,1), ji = 1, jpi,10) 
     538 
     539         WRITE(numout,*) 
     540         WRITE(numout,*) '          latitude and e2 scale factors' 
     541         WRITE(numout,*) '          -----------------------------' 
     542         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   & 
     543            &                     gphiv(1,jj), gphif(1,jj),   & 
     544            &                     e2t  (1,jj), e2u  (1,jj),   & 
     545            &                     e2v  (1,jj), jj = 1, jpj, 10 ) 
     546      ENDIF 
     547 
     548      !                                     ! ============================ 
     549      !                                     !        close the files  
     550      !                                     ! ============================ 
     551      CALL iom_close( inum ) 
     552      ! 
     5539300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    & 
     554            f19.10, 1x, f19.10, 1x, f19.10 ) 
     555   END SUBROUTINE dom_hgr 
     556 
     557 
    323558   SUBROUTINE dom_zgr 
    324559      !!---------------------------------------------------------------------- 
    325560      !!                ***  ROUTINE dom_zgr  *** 
    326561      !!                    
    327       !! ** Purpose :  set the depth of model levels and the resulting  
    328       !!      vertical scale factors. 
     562      !! ** Purpose :  Read the NetCDF file(s) which contain(s) all the 
     563      !!      ocean horizontal mesh informations and/or set the depth of model levels  
     564      !!      and the resulting vertical scale factors. 
    329565      !! 
    330566      !! ** Method  : - reference 1D vertical coordinate (gdep._1d, e3._1d) 
     
    338574      !! ** Action  :   define gdep., e3., mbathy and bathy 
    339575      !!---------------------------------------------------------------------- 
    340       INTEGER ::   ioptio = 0   ! temporary integer 
    341       INTEGER ::   ios 
     576      INTEGER  ::  ioptio = 0   ! temporary integer 
     577      INTEGER  ::  inum, ios 
     578      INTEGER  ::  ji, jj, jk, ik 
     579      REAL(wp) ::  zrefdep 
    342580      !! 
    343581      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav 
     582      REAL(wp), POINTER, DIMENSION(:,:) :: zprt, zprw 
    344583      !!---------------------------------------------------------------------- 
    345584 
     
    372611      IF ( ioptio == 33 )   CALL ctl_stop( ' isf cavity with off line module not yet done    ' ) 
    373612 
    374    END SUBROUTINE dom_zgr 
    375  
    376    SUBROUTINE dom_ctl 
    377       !!---------------------------------------------------------------------- 
    378       !!                     ***  ROUTINE dom_ctl  *** 
    379       !! 
    380       !! ** Purpose :   Domain control. 
    381       !! 
    382       !! ** Method  :   compute and print extrema of masked scale factors 
    383       !! 
    384       !! History : 
    385       !!   8.5  !  02-08  (G. Madec)    Original code 
    386       !!---------------------------------------------------------------------- 
    387       !! * Local declarations 
    388       INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2 
    389       INTEGER, DIMENSION(2) ::   iloc      !  
    390       REAL(wp) ::   ze1min, ze1max, ze2min, ze2max 
    391       !!---------------------------------------------------------------------- 
    392  
    393       ! Extrema of the scale factors 
    394  
    395       IF(lwp)WRITE(numout,*) 
    396       IF(lwp)WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors' 
    397       IF(lwp)WRITE(numout,*) '~~~~~~~' 
    398  
    399       IF (lk_mpp) THEN 
    400          CALL mpp_minloc( e1t(:,:), tmask(:,:,1), ze1min, iimi1,ijmi1 ) 
    401          CALL mpp_minloc( e2t(:,:), tmask(:,:,1), ze2min, iimi2,ijmi2 ) 
    402          CALL mpp_maxloc( e1t(:,:), tmask(:,:,1), ze1max, iima1,ijma1 ) 
    403          CALL mpp_maxloc( e2t(:,:), tmask(:,:,1), ze2max, iima2,ijma2 ) 
     613      IF(lwp) WRITE(numout,*) '          one file in "mesh_mask.nc" ' 
     614      CALL iom_open( 'mesh_zgr', inum ) 
     615 
     616      CALL iom_get( inum, jpdom_unknown, 'gdept_1d', gdept_1d ) ! depth 
     617      CALL iom_get( inum, jpdom_unknown, 'gdepw_1d', gdepw_1d ) 
     618      IF( ln_zco .OR. ln_zps ) THEN 
     619         CALL iom_get( inum, jpdom_unknown, 'e3t_1d'  , e3t_1d   )    ! reference scale factors 
     620         CALL iom_get( inum, jpdom_unknown, 'e3w_1d'  , e3w_1d   ) 
     621      ENDIF 
     622 
     623!!gm BUG in s-coordinate this does not work! 
     624      ! deepest/shallowest W level Above/Below ~10m 
     625      zrefdep = 10._wp - ( 0.1_wp * MINVAL(e3w_1d) )                 ! ref. depth with tolerance (10% of minimum layer thickness) 
     626      nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m 
     627      nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m 
     628!!gm end bug 
     629 
     630      IF(lwp) THEN 
     631         WRITE(numout,*) 
     632         WRITE(numout,*) '              Reference z-coordinate depth and scale factors:' 
     633         WRITE(numout, "(9x,' level   gdept    gdepw     e3t      e3w  ')" ) 
     634         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk ) 
     635      ENDIF 
     636 
     637      DO jk = 1, jpk 
     638         IF( e3w_1d  (jk) <= 0._wp .OR. e3t_1d  (jk) <= 0._wp )   CALL ctl_stop( ' e3w_1d or e3t_1d =< 0 ' ) 
     639         IF( gdepw_1d(jk) <  0._wp .OR. gdept_1d(jk) <  0._wp )   CALL ctl_stop( ' gdepw_1d or gdept_1d < 0 ' ) 
     640      END DO 
     641 
     642      IF( lk_vvl ) THEN 
     643          CALL iom_get( inum, jpdom_data, 'e3t_0', e3t_0(:,:,:) ) 
     644          CALL iom_get( inum, jpdom_data, 'e3u_0', e3u_0(:,:,:) ) 
     645          CALL iom_get( inum, jpdom_data, 'e3v_0', e3v_0(:,:,:) ) 
     646          CALL iom_get( inum, jpdom_data, 'e3w_0', e3w_0(:,:,:) ) 
     647          CALL iom_get( inum, jpdom_data, 'gdept_0', gdept_0(:,:,:) ) 
     648          CALL iom_get( inum, jpdom_data, 'gdepw_0', gdepw_0(:,:,:) ) 
     649          ht_0(:,:) = 0.0_wp                       ! Reference ocean depth at  T-points 
     650          DO jk = 1, jpk 
     651             ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk) 
     652          END DO 
    404653      ELSE 
    405          ze1min = MINVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )     
    406          ze2min = MINVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )     
    407          ze1max = MAXVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )     
    408          ze2max = MAXVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )     
    409  
    410          iloc  = MINLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 ) 
    411          iimi1 = iloc(1) + nimpp - 1 
    412          ijmi1 = iloc(2) + njmpp - 1 
    413          iloc  = MINLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 ) 
    414          iimi2 = iloc(1) + nimpp - 1 
    415          ijmi2 = iloc(2) + njmpp - 1 
    416          iloc  = MAXLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 ) 
    417          iima1 = iloc(1) + nimpp - 1 
    418          ijma1 = iloc(2) + njmpp - 1 
    419          iloc  = MAXLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 ) 
    420          iima2 = iloc(1) + nimpp - 1 
    421          ijma2 = iloc(2) + njmpp - 1 
    422       ENDIF 
    423  
    424       IF(lwp) THEN 
    425          WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1 
    426          WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1 
    427          WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2 
    428          WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2 
    429       ENDIF 
    430  
    431    END SUBROUTINE dom_ctl 
    432  
    433    SUBROUTINE dom_grd 
    434       !!---------------------------------------------------------------------- 
    435       !!                  ***  ROUTINE dom_grd  *** 
    436       !!                    
    437       !! ** Purpose :  Read the NetCDF file(s) which contain(s) all the 
    438       !!      ocean domain informations (mesh and mask arrays). This (these) 
    439       !!      file(s) is (are) used for visualisation (SAXO software) and 
    440       !!      diagnostic computation. 
    441       !! 
    442       !! ** Method  :   Read in a file all the arrays generated in routines 
    443       !!      domhgr, domzgr, and dommsk. Note: the file contain depends on 
    444       !!      the vertical coord. used (z-coord, partial steps, s-coord) 
    445       !!                    nmsh = 1  :   'mesh_mask.nc' file 
    446       !!                         = 2  :   'mesh.nc' and mask.nc' files 
    447       !!                         = 3  :   'mesh_hgr.nc', 'mesh_zgr.nc' and 
    448       !!                                  'mask.nc' files 
    449       !!      For huge size domain, use option 2 or 3 depending on your  
    450       !!      vertical coordinate. 
    451       !! 
    452       !! ** input file :  
    453       !!      meshmask.nc  : domain size, horizontal grid-point position, 
    454       !!                     masks, depth and vertical scale factors 
    455       !!---------------------------------------------------------------------- 
    456       USE iom 
    457       !! 
    458       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    459       INTEGER  ::   ik, inum0 , inum1 , inum2 , inum3 , inum4   ! local integers 
    460       REAL(wp) ::   zrefdep         ! local real 
    461       REAL(wp), POINTER, DIMENSION(:,:) :: zmbk, zprt, zprw 
    462       !!---------------------------------------------------------------------- 
    463  
    464       IF(lwp) WRITE(numout,*) 
    465       IF(lwp) WRITE(numout,*) 'dom_rea : read NetCDF mesh and mask information file(s)' 
    466       IF(lwp) WRITE(numout,*) '~~~~~~~' 
    467  
    468       CALL wrk_alloc( jpi, jpj, zmbk, zprt, zprw ) 
    469  
    470       zmbk(:,:) = 0._wp 
    471  
    472       SELECT CASE (nmsh) 
    473          !                                     ! ============================ 
    474          CASE ( 1 )                            !  create 'mesh_mask.nc' file 
    475             !                                  ! ============================ 
    476  
    477             IF(lwp) WRITE(numout,*) '          one file in "mesh_mask.nc" ' 
    478             CALL iom_open( 'mesh_mask', inum0 ) 
    479  
    480             inum2 = inum0                                            ! put all the informations 
    481             inum3 = inum0                                            ! in unit inum0 
    482             inum4 = inum0 
    483  
    484             !                                  ! ============================ 
    485          CASE ( 2 )                            !  create 'mesh.nc' and  
    486             !                                  !         'mask.nc' files 
    487             !                                  ! ============================ 
    488  
    489             IF(lwp) WRITE(numout,*) '          two files in "mesh.nc" and "mask.nc" ' 
    490             CALL iom_open( 'mesh', inum1 ) 
    491             CALL iom_open( 'mask', inum2 ) 
    492  
    493             inum3 = inum1                                            ! put mesh informations  
    494             inum4 = inum1                                            ! in unit inum1  
    495  
    496             !                                  ! ============================ 
    497          CASE ( 3 )                            !  create 'mesh_hgr.nc' 
    498             !                                  !         'mesh_zgr.nc' and 
    499             !                                  !         'mask.nc'     files 
    500             !                                  ! ============================ 
    501  
    502             IF(lwp) WRITE(numout,*) '          three files in "mesh_hgr.nc" , "mesh_zgr.nc" and "mask.nc" ' 
    503             CALL iom_open( 'mesh_hgr', inum3 ) ! create 'mesh_hgr.nc' 
    504             CALL iom_open( 'mesh_zgr', inum4 ) ! create 'mesh_zgr.nc' 
    505             CALL iom_open( 'mask'    , inum2 ) ! create 'mask.nc' 
    506  
    507             !                                  ! =========================== 
    508          CASE DEFAULT                          ! return error  
    509             !                                  ! mesh has to be provided 
    510             !                                  ! =========================== 
    511             CALL ctl_stop( ' OFFLINE mode requires the input mesh mask(s). ',   & 
    512             &                                 ' Invalid nn_msh value in the namelist (0 is not allowed)' ) 
    513  
    514          END SELECT 
    515  
    516          !                                                         ! masks (inum2)  
    517          CALL iom_get( inum2, jpdom_data, 'tmask', tmask ) 
    518          CALL iom_get( inum2, jpdom_data, 'umask', umask ) 
    519          CALL iom_get( inum2, jpdom_data, 'vmask', vmask ) 
    520          CALL iom_get( inum2, jpdom_data, 'fmask', fmask ) 
    521  
    522          CALL lbc_lnk( tmask, 'T', 1._wp )    ! Lateral boundary conditions 
    523          CALL lbc_lnk( umask, 'U', 1._wp )       
    524          CALL lbc_lnk( vmask, 'V', 1._wp ) 
    525          CALL lbc_lnk( fmask, 'F', 1._wp ) 
    526  
    527 #if defined key_c1d 
    528          ! set umask and vmask equal tmask in 1D configuration 
    529          IF(lwp) WRITE(numout,*) 
    530          IF(lwp) WRITE(numout,*) '**********  1D configuration : set umask and vmask equal tmask ********' 
    531          IF(lwp) WRITE(numout,*) '**********                                                     ********' 
    532  
    533          umask(:,:,:) = tmask(:,:,:) 
    534          vmask(:,:,:) = tmask(:,:,:) 
    535 #endif 
    536  
    537 #if defined key_degrad 
    538          CALL iom_get( inum2, jpdom_data, 'facvolt', facvol ) 
    539 #endif 
    540  
    541          !                                                         ! horizontal mesh (inum3) 
    542          CALL iom_get( inum3, jpdom_data, 'glamt', glamt ) 
    543          CALL iom_get( inum3, jpdom_data, 'glamu', glamu ) 
    544          CALL iom_get( inum3, jpdom_data, 'glamv', glamv ) 
    545          CALL iom_get( inum3, jpdom_data, 'glamf', glamf ) 
    546  
    547          CALL iom_get( inum3, jpdom_data, 'gphit', gphit ) 
    548          CALL iom_get( inum3, jpdom_data, 'gphiu', gphiu ) 
    549          CALL iom_get( inum3, jpdom_data, 'gphiv', gphiv ) 
    550          CALL iom_get( inum3, jpdom_data, 'gphif', gphif ) 
    551  
    552          CALL iom_get( inum3, jpdom_data, 'e1t', e1t ) 
    553          CALL iom_get( inum3, jpdom_data, 'e1u', e1u ) 
    554          CALL iom_get( inum3, jpdom_data, 'e1v', e1v ) 
    555           
    556          CALL iom_get( inum3, jpdom_data, 'e2t', e2t ) 
    557          CALL iom_get( inum3, jpdom_data, 'e2u', e2u ) 
    558          CALL iom_get( inum3, jpdom_data, 'e2v', e2v ) 
    559  
    560          CALL iom_get( inum3, jpdom_data, 'ff', ff ) 
    561  
    562          CALL iom_get( inum4, jpdom_data, 'mbathy', zmbk )              ! number of ocean t-points 
    563          mbathy (:,:) = INT( zmbk(:,:) ) 
    564          misfdep(:,:) = 1                                               ! ice shelf case not yet done 
    565           
    566          CALL zgr_bot_level                                             ! mbk. arrays (deepest ocean t-, u- & v-points 
    567          ! 
    568654         IF( ln_sco ) THEN                                         ! s-coordinate 
    569             CALL iom_get( inum4, jpdom_data, 'hbatt', hbatt ) 
    570             CALL iom_get( inum4, jpdom_data, 'hbatu', hbatu ) 
    571             CALL iom_get( inum4, jpdom_data, 'hbatv', hbatv ) 
    572             CALL iom_get( inum4, jpdom_data, 'hbatf', hbatf ) 
     655            CALL iom_get( inum, jpdom_data, 'hbatt', hbatt ) 
     656            CALL iom_get( inum, jpdom_data, 'hbatu', hbatu ) 
     657            CALL iom_get( inum, jpdom_data, 'hbatv', hbatv ) 
     658            CALL iom_get( inum, jpdom_data, 'hbatf', hbatf ) 
    573659             
    574             CALL iom_get( inum4, jpdom_unknown, 'gsigt', gsigt ) ! scaling coef. 
    575             CALL iom_get( inum4, jpdom_unknown, 'gsigw', gsigw ) 
    576             CALL iom_get( inum4, jpdom_unknown, 'gsi3w', gsi3w )  
    577             CALL iom_get( inum4, jpdom_unknown, 'esigt', esigt ) 
    578             CALL iom_get( inum4, jpdom_unknown, 'esigw', esigw ) 
    579  
    580             CALL iom_get( inum4, jpdom_data, 'e3t_0', fse3t_n(:,:,:) ) ! scale factors 
    581             CALL iom_get( inum4, jpdom_data, 'e3u_0', fse3u_n(:,:,:) ) 
    582             CALL iom_get( inum4, jpdom_data, 'e3v_0', fse3v_n(:,:,:) ) 
    583             CALL iom_get( inum4, jpdom_data, 'e3w_0', fse3w_n(:,:,:) ) 
    584  
    585             CALL iom_get( inum4, jpdom_unknown, 'gdept_1d', gdept_1d ) ! depth 
    586             CALL iom_get( inum4, jpdom_unknown, 'gdepw_1d', gdepw_1d ) 
     660            CALL iom_get( inum, jpdom_unknown, 'gsigt', gsigt ) ! scaling coef. 
     661            CALL iom_get( inum, jpdom_unknown, 'gsigw', gsigw ) 
     662            CALL iom_get( inum, jpdom_unknown, 'gsi3w', gsi3w )  
     663            CALL iom_get( inum, jpdom_unknown, 'esigt', esigt ) 
     664            CALL iom_get( inum, jpdom_unknown, 'esigw', esigw ) 
     665 
     666            CALL iom_get( inum, jpdom_data, 'e3t_0', fse3t_n(:,:,:) ) ! scale factors 
     667            CALL iom_get( inum, jpdom_data, 'e3u_0', fse3u_n(:,:,:) ) 
     668            CALL iom_get( inum, jpdom_data, 'e3v_0', fse3v_n(:,:,:) ) 
     669            CALL iom_get( inum, jpdom_data, 'e3w_0', fse3w_n(:,:,:) ) 
    587670         ENDIF 
    588671 
    589672  
    590673         IF( ln_zps ) THEN                                           ! z-coordinate - partial steps 
    591             CALL iom_get( inum4, jpdom_unknown, 'gdept_1d', gdept_1d )  ! reference depth 
    592             CALL iom_get( inum4, jpdom_unknown, 'gdepw_1d', gdepw_1d ) 
    593             CALL iom_get( inum4, jpdom_unknown, 'e3t_1d'  , e3t_1d   )    ! reference scale factors 
    594             CALL iom_get( inum4, jpdom_unknown, 'e3w_1d'  , e3w_1d   ) 
    595674            ! 
    596             IF( nmsh <= 6 ) THEN                                        ! 3D vertical scale factors 
    597                CALL iom_get( inum4, jpdom_data, 'e3t_0', fse3t_n(:,:,:) ) 
    598                CALL iom_get( inum4, jpdom_data, 'e3u_0', fse3u_n(:,:,:) ) 
    599                CALL iom_get( inum4, jpdom_data, 'e3v_0', fse3v_n(:,:,:) ) 
    600                CALL iom_get( inum4, jpdom_data, 'e3w_0', fse3w_n(:,:,:) ) 
     675            IF( iom_varid( inum, 'e3t_0', ldstop = .FALSE. ) > 0 ) THEN 
     676               CALL iom_get( inum, jpdom_data, 'e3t_0', fse3t_n(:,:,:) ) 
     677               CALL iom_get( inum, jpdom_data, 'e3u_0', fse3u_n(:,:,:) ) 
     678               CALL iom_get( inum, jpdom_data, 'e3v_0', fse3v_n(:,:,:) ) 
     679               CALL iom_get( inum, jpdom_data, 'e3w_0', fse3w_n(:,:,:) ) 
    601680            ELSE                                                        ! 2D bottom scale factors 
    602                CALL iom_get( inum4, jpdom_data, 'e3t_ps', e3tp ) 
    603                CALL iom_get( inum4, jpdom_data, 'e3w_ps', e3wp ) 
     681               CALL iom_get( inum, jpdom_data, 'e3t_ps', e3tp ) 
     682               CALL iom_get( inum, jpdom_data, 'e3w_ps', e3wp ) 
    604683               !                                                        ! deduces the 3D scale factors 
    605684               DO jk = 1, jpk 
     
    633712            END IF 
    634713 
    635             IF( iom_varid( inum4, 'gdept_0', ldstop = .FALSE. ) > 0 ) THEN   ! 3D depth of t- and w-level 
    636                CALL iom_get( inum4, jpdom_data, 'gdept_0', fsdept_n(:,:,:) ) 
    637                CALL iom_get( inum4, jpdom_data, 'gdepw_0', fsdepw_n(:,:,:) ) 
     714            IF( iom_varid( inum, 'gdept_0', ldstop = .FALSE. ) > 0 ) THEN   ! 3D depth of t- and w-level 
     715               CALL iom_get( inum, jpdom_data, 'gdept_0', fsdept_n(:,:,:) ) 
     716               CALL iom_get( inum, jpdom_data, 'gdepw_0', fsdepw_n(:,:,:) ) 
    638717            ELSE                                                           ! 2D bottom depth 
    639                CALL iom_get( inum4, jpdom_data, 'hdept', zprt ) 
    640                CALL iom_get( inum4, jpdom_data, 'hdepw', zprw ) 
     718               CALL wrk_alloc( jpi, jpj, zprt, zprw ) 
     719               ! 
     720               CALL iom_get( inum, jpdom_data, 'hdept', zprt ) 
     721               CALL iom_get( inum, jpdom_data, 'hdepw', zprw ) 
    641722               ! 
    642723               DO jk = 1, jpk                                              ! deduces the 3D depth 
     
    654735                  END DO 
    655736               END DO 
     737               CALL wrk_dealloc( jpi, jpj, zprt, zprw ) 
    656738            ENDIF 
    657739            ! 
     
    659741 
    660742         IF( ln_zco ) THEN           ! Vertical coordinates and scales factors 
    661             CALL iom_get( inum4, jpdom_unknown, 'gdept_1d', gdept_1d ) ! depth 
    662             CALL iom_get( inum4, jpdom_unknown, 'gdepw_1d', gdepw_1d ) 
    663             CALL iom_get( inum4, jpdom_unknown, 'e3t_1d'  , e3t_1d   ) 
    664             CALL iom_get( inum4, jpdom_unknown, 'e3w_1d'  , e3w_1d   ) 
    665743            DO jk = 1, jpk 
    666744               fse3t_n(:,:,jk) = e3t_1d(jk)                              ! set to the ref. factors 
     
    672750            END DO 
    673751         ENDIF 
    674  
    675 !!gm BUG in s-coordinate this does not work! 
    676       ! deepest/shallowest W level Above/Below ~10m 
    677       zrefdep = 10._wp - ( 0.1_wp * MINVAL(e3w_1d) )                 ! ref. depth with tolerance (10% of minimum layer thickness) 
    678       nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m 
    679       nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m 
    680 !!gm end bug 
    681  
    682       ! Control printing : Grid informations (if not restart) 
    683       ! ---------------- 
    684  
    685       IF(lwp .AND. .NOT.ln_rstart ) THEN 
    686          WRITE(numout,*) 
    687          WRITE(numout,*) '          longitude and e1 scale factors' 
    688          WRITE(numout,*) '          ------------------------------' 
    689          WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   & 
    690             glamv(ji,1), glamf(ji,1),   & 
    691             e1t(ji,1), e1u(ji,1),   & 
    692             e1v(ji,1), ji = 1, jpi,10) 
    693 9300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    & 
    694             f19.10, 1x, f19.10, 1x, f19.10 ) 
    695  
    696          WRITE(numout,*) 
    697          WRITE(numout,*) '          latitude and e2 scale factors' 
    698          WRITE(numout,*) '          -----------------------------' 
    699          WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   & 
    700             &                     gphiv(1,jj), gphif(1,jj),   & 
    701             &                     e2t  (1,jj), e2u  (1,jj),   & 
    702             &                     e2v  (1,jj), jj = 1, jpj, 10 ) 
    703       ENDIF 
    704  
    705  
    706       IF( nprint == 1 .AND. lwp ) THEN 
    707          WRITE(numout,*) '          e1u e2u ' 
    708          CALL prihre( e1u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout ) 
    709          CALL prihre( e2u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout ) 
    710          WRITE(numout,*) '          e1v e2v  ' 
    711          CALL prihre( e1v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout ) 
    712          CALL prihre( e2v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout ) 
    713       ENDIF 
    714  
    715       IF(lwp) THEN 
    716          WRITE(numout,*) 
    717          WRITE(numout,*) '              Reference z-coordinate depth and scale factors:' 
    718          WRITE(numout, "(9x,' level   gdept    gdepw     e3t      e3w  ')" ) 
    719          WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk ) 
    720       ENDIF 
    721  
    722       DO jk = 1, jpk 
    723          IF( e3w_1d  (jk) <= 0._wp .OR. e3t_1d  (jk) <= 0._wp )   CALL ctl_stop( ' e3w_1d or e3t_1d =< 0 ' ) 
    724          IF( gdepw_1d(jk) <  0._wp .OR. gdept_1d(jk) <  0._wp )   CALL ctl_stop( ' gdepw_1d or gdept_1d < 0 ' ) 
    725       END DO 
     752         ! 
     753      ENDIF 
    726754      !                                     ! ============================ 
    727755      !                                     !        close the files  
    728756      !                                     ! ============================ 
    729       SELECT CASE ( nmsh ) 
    730          CASE ( 1 )                 
    731             CALL iom_close( inum0 ) 
    732          CASE ( 2 ) 
    733             CALL iom_close( inum1 ) 
    734             CALL iom_close( inum2 ) 
    735          CASE ( 3 ) 
    736             CALL iom_close( inum2 ) 
    737             CALL iom_close( inum3 ) 
    738             CALL iom_close( inum4 ) 
    739       END SELECT 
    740       ! 
    741       CALL wrk_dealloc( jpi, jpj, zmbk, zprt, zprw ) 
    742       ! 
    743    END SUBROUTINE dom_grd 
    744  
    745  
    746    SUBROUTINE zgr_bot_level 
    747       !!---------------------------------------------------------------------- 
    748       !!                    ***  ROUTINE zgr_bot_level  *** 
    749       !! 
    750       !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays) 
    751       !! 
    752       !! ** Method  :   computes from mbathy with a minimum value of 1 over land 
    753       !! 
    754       !! ** Action  :   mbkt, mbku, mbkv :   vertical indices of the deeptest 
    755       !!                                     ocean level at t-, u- & v-points 
    756       !!                                     (min value = 1 over land) 
    757       !!---------------------------------------------------------------------- 
    758       ! 
    759       INTEGER ::   ji, jj   ! dummy loop indices 
    760       REAL(wp), POINTER, DIMENSION(:,:) :: zmbk 
    761       !!---------------------------------------------------------------------- 
    762  
    763       ! 
    764       IF(lwp) WRITE(numout,*) 
    765       IF(lwp) WRITE(numout,*) '    zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels ' 
    766       IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~' 
    767       ! 
    768       CALL wrk_alloc( jpi, jpj, zmbk ) 
    769       ! 
    770       mbkt(:,:) = MAX( mbathy(:,:) , 1 )    ! bottom k-index of T-level (=1 over land) 
    771       mikt(:,:) = 1 ; miku(:,:) = 1; mikv(:,:) = 1; ! top k-index of T-level (=1 over open ocean; >1 beneath ice shelf) 
    772       !                                     ! bottom k-index of W-level = mbkt+1 
    773       DO jj = 1, jpjm1                      ! bottom k-index of u- (v-) level 
    774          DO ji = 1, jpim1 
    775             mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  ) 
    776             mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  ) 
    777          END DO 
    778       END DO 
    779       ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk  
    780       zmbk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 
    781       zmbk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 
    782       ! 
    783       CALL wrk_dealloc( jpi, jpj, zmbk ) 
    784       ! 
    785    END SUBROUTINE zgr_bot_level 
    786  
    787    SUBROUTINE dom_msk 
    788       !!--------------------------------------------------------------------- 
    789       !!                 ***  ROUTINE dom_msk  *** 
    790       !! 
    791       !! ** Purpose :   Off-line case: defines the interior domain T-mask. 
    792       !! 
    793       !! ** Method  :   The interior ocean/land mask is computed from tmask 
    794       !!              setting to zero the duplicated row and lines due to 
    795       !!              MPP exchange halos, est-west cyclic and north fold 
    796       !!              boundary conditions. 
    797       !! 
    798       !! ** Action :   tmask_i  : interiorland/ocean mask at t-point 
    799       !!               tpol     : ??? 
    800       !!---------------------------------------------------------------------- 
    801       ! 
    802       INTEGER  ::   ji, jj, jk                   ! dummy loop indices 
    803       INTEGER  ::   iif, iil, ijf, ijl       ! local integers 
    804       INTEGER, POINTER, DIMENSION(:,:) ::  imsk  
    805       ! 
    806       !!--------------------------------------------------------------------- 
    807        
    808       CALL wrk_alloc( jpi, jpj, imsk ) 
    809       ! 
    810       ! Interior domain mask (used for global sum) 
    811       ! -------------------- 
    812       ssmask(:,:)  = tmask(:,:,1) 
    813       tmask_i(:,:) = tmask(:,:,1) 
    814       iif = jpreci                        ! thickness of exchange halos in i-axis 
    815       iil = nlci - jpreci + 1 
    816       ijf = jprecj                        ! thickness of exchange halos in j-axis 
    817       ijl = nlcj - jprecj + 1 
    818       ! 
    819       tmask_i( 1 :iif,   :   ) = 0._wp    ! first columns 
    820       tmask_i(iil:jpi,   :   ) = 0._wp    ! last  columns (including mpp extra columns) 
    821       tmask_i(   :   , 1 :ijf) = 0._wp    ! first rows 
    822       tmask_i(   :   ,ijl:jpj) = 0._wp    ! last  rows (including mpp extra rows) 
    823       ! 
    824       !                                   ! north fold mask 
    825       tpol(1:jpiglo) = 1._wp 
    826       !                                 
    827       IF( jperio == 3 .OR. jperio == 4 )   tpol(jpiglo/2+1:jpiglo) = 0._wp    ! T-point pivot 
    828       IF( jperio == 5 .OR. jperio == 6 )   tpol(     1    :jpiglo) = 0._wp    ! F-point pivot 
    829       IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot: only half of the nlcj-1 row 
    830          IF( mjg(ijl-1) == jpjglo-1 ) THEN 
    831             DO ji = iif+1, iil-1 
    832                tmask_i(ji,ijl-1) = tmask_i(ji,ijl-1) * tpol(mig(ji)) 
    833             END DO 
    834          ENDIF 
    835       ENDIF  
    836       ! 
    837       ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at 
    838       ! least 1 wet u point 
    839       DO jj = 1, jpjm1 
    840          DO ji = 1, fs_jpim1   ! vector loop 
    841             umask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:))) 
    842             vmask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:))) 
    843          END DO 
    844          DO ji = 1, jpim1      ! NO vector opt. 
    845             fmask_i(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   & 
    846                &            * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:))) 
    847          END DO 
    848       END DO 
    849       CALL lbc_lnk( umask_i, 'U', 1._wp )      ! Lateral boundary conditions 
    850       CALL lbc_lnk( vmask_i, 'V', 1._wp ) 
    851       CALL lbc_lnk( fmask_i, 'F', 1._wp ) 
    852  
    853       ! 3. Ocean/land mask at wu-, wv- and w points  
    854       !---------------------------------------------- 
    855       wmask (:,:,1) = tmask(:,:,1) ! ???????? 
    856       wumask(:,:,1) = umask(:,:,1) ! ???????? 
    857       wvmask(:,:,1) = vmask(:,:,1) ! ???????? 
    858       DO jk=2,jpk 
    859          wmask (:,:,jk)=tmask(:,:,jk) * tmask(:,:,jk-1) 
    860          wumask(:,:,jk)=umask(:,:,jk) * umask(:,:,jk-1)    
    861          wvmask(:,:,jk)=vmask(:,:,jk) * vmask(:,:,jk-1) 
    862       END DO 
    863       ! 
    864       IF( nprint == 1 .AND. lwp ) THEN    ! Control print 
    865          imsk(:,:) = INT( tmask_i(:,:) ) 
    866          WRITE(numout,*) ' tmask_i : ' 
    867          CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 1, numout) 
    868          WRITE (numout,*) 
    869          WRITE (numout,*) ' dommsk: tmask for each level' 
    870          WRITE (numout,*) ' ----------------------------' 
    871          DO jk = 1, jpk 
    872             imsk(:,:) = INT( tmask(:,:,jk) ) 
    873             WRITE(numout,*) 
    874             WRITE(numout,*) ' level = ',jk 
    875             CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 1, numout) 
    876          END DO 
    877       ENDIF 
    878       ! 
    879       CALL wrk_dealloc( jpi, jpj, imsk ) 
    880       ! 
    881    END SUBROUTINE dom_msk 
     757      CALL iom_close( inum ) 
     758      ! 
     759      ! 
     760   END SUBROUTINE dom_zgr 
     761 
     762   SUBROUTINE dom_ctl 
     763      !!---------------------------------------------------------------------- 
     764      !!                     ***  ROUTINE dom_ctl  *** 
     765      !! 
     766      !! ** Purpose :   Domain control. 
     767      !! 
     768      !! ** Method  :   compute and print extrema of masked scale factors 
     769      !! 
     770      !! History : 
     771      !!   8.5  !  02-08  (G. Madec)    Original code 
     772      !!---------------------------------------------------------------------- 
     773      !! * Local declarations 
     774      INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2 
     775      INTEGER, DIMENSION(2) ::   iloc      !  
     776      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max 
     777      !!---------------------------------------------------------------------- 
     778 
     779      ! Extrema of the scale factors 
     780 
     781      IF(lwp)WRITE(numout,*) 
     782      IF(lwp)WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors' 
     783      IF(lwp)WRITE(numout,*) '~~~~~~~' 
     784 
     785      IF (lk_mpp) THEN 
     786         CALL mpp_minloc( e1t(:,:), tmask(:,:,1), ze1min, iimi1,ijmi1 ) 
     787         CALL mpp_minloc( e2t(:,:), tmask(:,:,1), ze2min, iimi2,ijmi2 ) 
     788         CALL mpp_maxloc( e1t(:,:), tmask(:,:,1), ze1max, iima1,ijma1 ) 
     789         CALL mpp_maxloc( e2t(:,:), tmask(:,:,1), ze2max, iima2,ijma2 ) 
     790      ELSE 
     791         ze1min = MINVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )     
     792         ze2min = MINVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )     
     793         ze1max = MAXVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )     
     794         ze2max = MAXVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )     
     795 
     796         iloc  = MINLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 ) 
     797         iimi1 = iloc(1) + nimpp - 1 
     798         ijmi1 = iloc(2) + njmpp - 1 
     799         iloc  = MINLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 ) 
     800         iimi2 = iloc(1) + nimpp - 1 
     801         ijmi2 = iloc(2) + njmpp - 1 
     802         iloc  = MAXLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 ) 
     803         iima1 = iloc(1) + nimpp - 1 
     804         ijma1 = iloc(2) + njmpp - 1 
     805         iloc  = MAXLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 ) 
     806         iima2 = iloc(1) + nimpp - 1 
     807         ijma2 = iloc(2) + njmpp - 1 
     808      ENDIF 
     809 
     810      IF(lwp) THEN 
     811         WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1 
     812         WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1 
     813         WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2 
     814         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2 
     815      ENDIF 
     816 
     817   END SUBROUTINE dom_ctl 
    882818 
    883819   !!====================================================================== 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90

    r5767 r7522  
    2222   USE c1d             ! 1D configuration: lk_c1d 
    2323   USE dom_oce         ! ocean domain: variables 
     24   USE domvvl          ! variable volume 
    2425   USE zdf_oce         ! ocean vertical physics: variables 
    2526   USE sbc_oce         ! surface module: variables 
     
    2829   USE trabbl          ! active tracer: bottom boundary layer 
    2930   USE ldfslp          ! lateral diffusion: iso-neutral slopes 
     31   USE sbcrnf          ! river runoffs 
    3032   USE ldfeiv          ! eddy induced velocity coef.  
    3133   USE ldftra_oce      ! ocean tracer   lateral physics 
     
    3941   USE prtctl          ! print control 
    4042   USE fldread         ! read input fields  
     43   USE wrk_nemo        ! Memory allocation  
    4144   USE timing          ! Timing 
     45   USE trc, ONLY : ln_rsttr, numrtr, numrtw, lrst_trc 
    4246 
    4347   IMPLICIT NONE 
     
    4650   PUBLIC   dta_dyn_init   ! called by opa.F90 
    4751   PUBLIC   dta_dyn        ! called by step.F90 
    48  
    49    CHARACTER(len=100) ::   cn_dir       !: Root directory for location of ssr files 
    50    LOGICAL            ::   ln_dynwzv    !: vertical velocity read in a file (T) or computed from u/v (F) 
    51    LOGICAL            ::   ln_dynbbl    !: bbl coef read in a file (T) or computed (F) 
    52    LOGICAL            ::   ln_degrad    !: degradation option enabled or not 
    53    LOGICAL            ::   ln_dynrnf    !: read runoff data in file (T) or set to zero (F) 
    54  
    55    INTEGER  , PARAMETER ::   jpfld = 21     ! maximum number of fields to read 
     52   PUBLIC   dta_dyn_swp   ! called by step.F90 
     53 
     54   CHARACTER(len=100) ::   cn_dir          !: Root directory for location of ssr files 
     55   LOGICAL            ::   ln_dynrnf       !: read runoff data in file (T) or set to zero (F) 
     56   LOGICAL            ::   ln_dynrnf_depth       !: read runoff data in file (T) or set to zero (F) 
     57   REAL(wp)           ::   fwbcorr 
     58 
     59 
     60   INTEGER  , PARAMETER ::   jpfld = 20     ! maximum number of fields to read 
    5661   INTEGER  , SAVE      ::   jf_tem         ! index of temperature 
    5762   INTEGER  , SAVE      ::   jf_sal         ! index of salinity 
    58    INTEGER  , SAVE      ::   jf_uwd         ! index of u-wind 
    59    INTEGER  , SAVE      ::   jf_vwd         ! index of v-wind 
    60    INTEGER  , SAVE      ::   jf_wwd         ! index of w-wind 
     63   INTEGER  , SAVE      ::   jf_uwd         ! index of u-transport 
     64   INTEGER  , SAVE      ::   jf_vwd         ! index of v-transport 
     65   INTEGER  , SAVE      ::   jf_wwd         ! index of v-transport 
    6166   INTEGER  , SAVE      ::   jf_avt         ! index of Kz 
    6267   INTEGER  , SAVE      ::   jf_mld         ! index of mixed layer deptht 
    6368   INTEGER  , SAVE      ::   jf_emp         ! index of water flux 
     69   INTEGER  , SAVE      ::   jf_empb        ! index of water flux 
    6470   INTEGER  , SAVE      ::   jf_qsr         ! index of solar radiation 
    6571   INTEGER  , SAVE      ::   jf_wnd         ! index of wind speed 
    6672   INTEGER  , SAVE      ::   jf_ice         ! index of sea ice cover 
    6773   INTEGER  , SAVE      ::   jf_rnf         ! index of river runoff 
     74   INTEGER  , SAVE      ::   jf_fmf         ! index of downward salt flux 
    6875   INTEGER  , SAVE      ::   jf_ubl         ! index of u-bbl coef 
    6976   INTEGER  , SAVE      ::   jf_vbl         ! index of v-bbl coef 
    70    INTEGER  , SAVE      ::   jf_ahu         ! index of u-diffusivity coef 
    71    INTEGER  , SAVE      ::   jf_ahv         ! index of v-diffusivity coef  
    72    INTEGER  , SAVE      ::   jf_ahw         ! index of w-diffusivity coef 
    73    INTEGER  , SAVE      ::   jf_eiu         ! index of u-eiv 
    74    INTEGER  , SAVE      ::   jf_eiv         ! index of v-eiv 
    75    INTEGER  , SAVE      ::   jf_eiw         ! index of w-eiv 
    76    INTEGER  , SAVE      ::   jf_fmf         ! index of downward salt flux 
    77  
    78    TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_dyn  ! structure of input fields (file informations, fields read) 
     77   INTEGER  , SAVE      ::   jf_div         ! index of e3t 
     78 
     79 
     80   TYPE(FLD), ALLOCATABLE, SAVE, DIMENSION(:) :: sf_dyn  ! structure of input fields (file informations, fields read) 
    7981   !                                               !  
    80    REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wdta       ! vertical velocity at 2 time step 
    81    REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:  ) :: wnow       ! vertical velocity at 2 time step 
    8282   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: uslpdta    ! zonal isopycnal slopes 
    8383   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: vslpdta    ! meridional isopycnal slopes 
    8484   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpidta   ! zonal diapycnal slopes 
    8585   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpjdta   ! meridional diapycnal slopes 
    86    REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: uslpnow    ! zonal isopycnal slopes 
    87    REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: vslpnow    ! meridional isopycnal slopes 
    88    REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: wslpinow   ! zonal diapycnal slopes 
    89    REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: wslpjnow   ! meridional diapycnal slopes 
    90  
    91    INTEGER :: nrecprev_tem , nrecprev_uwd 
     86 
     87   INTEGER, SAVE  :: nprevrec, nsecdyn 
    9288 
    9389   !! * Substitutions 
     
    113109      !!---------------------------------------------------------------------- 
    114110      ! 
    115       USE oce, ONLY:  zts    => tsa  
    116       USE oce, ONLY:  zuslp  => ua   , zvslp  => va 
    117       USE oce, ONLY:  zwslpi => rotb , zwslpj => rotn 
    118       USE oce, ONLY:  zu     => ub   , zv     => vb,  zw => hdivb 
    119       ! 
     111      USE oce, ONLY:  zhdivtr => ua 
    120112      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    121       ! 
    122       INTEGER  ::   ji, jj     ! dummy loop indices 
    123       INTEGER  ::   isecsbc    ! number of seconds between Jan. 1st 00h of nit000 year and the middle of time step 
    124       REAL(wp) ::   ztinta     ! ratio applied to after  records when doing time interpolation 
    125       REAL(wp) ::   ztintb     ! ratio applied to before records when doing time interpolation 
    126       INTEGER  ::   iswap_tem, iswap_uwd    !  
     113      INTEGER             ::   ji, jj, jk 
     114      REAL(wp), POINTER, DIMENSION(:,:)   :: zemp 
     115      ! 
    127116      !!---------------------------------------------------------------------- 
    128117       
     
    130119      IF( nn_timing == 1 )  CALL timing_start( 'dta_dyn') 
    131120      ! 
    132       isecsbc = nsec_year + nsec1jan000  
    133       ! 
    134       IF( kt == nit000 ) THEN 
    135          nrecprev_tem = 0 
    136          nrecprev_uwd = 0 
    137          ! 
    138          CALL fld_read( kt, 1, sf_dyn )      !==   read data at kt time step   ==! 
    139          ! 
    140          IF( lk_ldfslp .AND. .NOT.lk_c1d .AND. sf_dyn(jf_tem)%ln_tint ) THEN    ! Computes slopes (here avt is used as workspace)                        
    141             zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,1) * tmask(:,:,:)   ! temperature 
    142             zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:)   ! salinity  
    143             avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,1) * tmask(:,:,:)   ! vertical diffusive coef. 
    144             CALL dta_dyn_slp( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
    145             uslpdta (:,:,:,1) = zuslp (:,:,:)  
    146             vslpdta (:,:,:,1) = zvslp (:,:,:)  
    147             wslpidta(:,:,:,1) = zwslpi(:,:,:)  
    148             wslpjdta(:,:,:,1) = zwslpj(:,:,:)  
    149          ENDIF 
    150          IF( ln_dynwzv .AND. sf_dyn(jf_uwd)%ln_tint )  THEN    ! compute vertical velocity from u/v 
    151             zu(:,:,:) = sf_dyn(jf_uwd)%fdta(:,:,:,1) 
    152             zv(:,:,:) = sf_dyn(jf_vwd)%fdta(:,:,:,1) 
    153             CALL dta_dyn_wzv( zu, zv, zw ) 
    154             wdta(:,:,:,1) = zw(:,:,:) * tmask(:,:,:) 
    155          ENDIF 
    156       ELSE 
    157          nrecprev_tem = sf_dyn(jf_tem)%nrec_a(2) 
    158          nrecprev_uwd = sf_dyn(jf_uwd)%nrec_a(2) 
    159          ! 
    160          CALL fld_read( kt, 1, sf_dyn )      !==   read data at kt time step   ==! 
    161          ! 
    162       ENDIF 
    163       !  
    164       IF( lk_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace)                        
    165          iswap_tem = 0 
    166          IF(  kt /= nit000 .AND. ( sf_dyn(jf_tem)%nrec_a(2) - nrecprev_tem ) /= 0 )  iswap_tem = 1 
    167          IF( ( isecsbc > sf_dyn(jf_tem)%nrec_b(2) .AND. iswap_tem == 1 ) .OR. kt == nit000 )  THEN    ! read/update the after data 
    168             IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt 
    169             IF( sf_dyn(jf_tem)%ln_tint ) THEN                 ! time interpolation of data 
    170                IF( kt /= nit000 ) THEN 
    171                   uslpdta (:,:,:,1) =  uslpdta (:,:,:,2)         ! swap the data 
    172                   vslpdta (:,:,:,1) =  vslpdta (:,:,:,2)   
    173                   wslpidta(:,:,:,1) =  wslpidta(:,:,:,2)  
    174                   wslpjdta(:,:,:,1) =  wslpjdta(:,:,:,2)  
    175                ENDIF 
    176                ! 
    177                zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:)   ! temperature 
    178                zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity  
    179                avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef. 
    180                CALL dta_dyn_slp( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
    181                ! 
    182                uslpdta (:,:,:,2) = zuslp (:,:,:)  
    183                vslpdta (:,:,:,2) = zvslp (:,:,:)  
    184                wslpidta(:,:,:,2) = zwslpi(:,:,:)  
    185                wslpjdta(:,:,:,2) = zwslpj(:,:,:)  
    186             ELSE 
    187                zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) 
    188                zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) 
    189                avt(:,:,:)        = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:) 
    190                CALL dta_dyn_slp( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
    191                uslpnow (:,:,:)   = zuslp (:,:,:)  
    192                vslpnow (:,:,:)   = zvslp (:,:,:)  
    193                wslpinow(:,:,:)   = zwslpi(:,:,:)  
    194                wslpjnow(:,:,:)   = zwslpj(:,:,:)  
    195             ENDIF 
    196          ENDIF 
    197          IF( sf_dyn(jf_tem)%ln_tint )  THEN 
    198             ztinta =  REAL( isecsbc - sf_dyn(jf_tem)%nrec_b(2), wp )  & 
    199                &    / REAL( sf_dyn(jf_tem)%nrec_a(2) - sf_dyn(jf_tem)%nrec_b(2), wp ) 
    200             ztintb =  1. - ztinta 
    201             uslp (:,:,:) = ztintb * uslpdta (:,:,:,1)  + ztinta * uslpdta (:,:,:,2)   
    202             vslp (:,:,:) = ztintb * vslpdta (:,:,:,1)  + ztinta * vslpdta (:,:,:,2)   
    203             wslpi(:,:,:) = ztintb * wslpidta(:,:,:,1)  + ztinta * wslpidta(:,:,:,2)   
    204             wslpj(:,:,:) = ztintb * wslpjdta(:,:,:,1)  + ztinta * wslpjdta(:,:,:,2)   
    205          ELSE 
    206             uslp (:,:,:) = uslpnow (:,:,:) 
    207             vslp (:,:,:) = vslpnow (:,:,:) 
    208             wslpi(:,:,:) = wslpinow(:,:,:) 
    209             wslpj(:,:,:) = wslpjnow(:,:,:) 
    210          ENDIF 
    211       ENDIF 
    212       ! 
    213       IF( ln_dynwzv )  THEN    ! compute vertical velocity from u/v 
    214          iswap_uwd = 0 
    215          IF(  kt /= nit000 .AND. ( sf_dyn(jf_uwd)%nrec_a(2) - nrecprev_uwd ) /= 0 )  iswap_uwd = 1 
    216          IF( ( isecsbc > sf_dyn(jf_uwd)%nrec_b(2) .AND. iswap_uwd == 1 ) .OR. kt == nit000 )  THEN    ! read/update the after data 
    217             IF(lwp) WRITE(numout,*) ' Compute new vertical velocity at kt = ', kt 
    218             IF(lwp) WRITE(numout,*) 
    219             IF( sf_dyn(jf_uwd)%ln_tint ) THEN                 ! time interpolation of data 
    220                IF( kt /= nit000 )  THEN 
    221                   wdta(:,:,:,1) =  wdta(:,:,:,2)     ! swap the data for initialisation 
    222                ENDIF 
    223                zu(:,:,:) = sf_dyn(jf_uwd)%fdta(:,:,:,2) 
    224                zv(:,:,:) = sf_dyn(jf_vwd)%fdta(:,:,:,2) 
    225                CALL dta_dyn_wzv( zu, zv, zw ) 
    226                wdta(:,:,:,2) = zw(:,:,:) * tmask(:,:,:) 
    227             ELSE 
    228                zu(:,:,:) = sf_dyn(jf_uwd)%fnow(:,:,:)  
    229                zv(:,:,:) = sf_dyn(jf_vwd)%fnow(:,:,:) 
    230                CALL dta_dyn_wzv( zu, zv, zw ) 
    231                wnow(:,:,:)  = zw(:,:,:) * tmask(:,:,:) 
    232             ENDIF 
    233          ENDIF 
    234          IF( sf_dyn(jf_uwd)%ln_tint )  THEN 
    235             ztinta =  REAL( isecsbc - sf_dyn(jf_uwd)%nrec_b(2), wp )  & 
    236                &    / REAL( sf_dyn(jf_uwd)%nrec_a(2) - sf_dyn(jf_uwd)%nrec_b(2), wp ) 
    237             ztintb =  1. - ztinta 
    238             wn(:,:,:) = ztintb * wdta(:,:,:,1)  + ztinta * wdta(:,:,:,2)   
    239          ELSE 
    240             wn(:,:,:) = wnow(:,:,:) 
    241          ENDIF 
    242       ENDIF 
     121      ! 
     122      nsecdyn = nsec_year + nsec1jan000   ! number of seconds between Jan. 1st 00h of nit000 year and the middle of time step 
     123      ! 
     124      IF( kt == nit000 ) THEN    ;    nprevrec = 0 
     125      ELSE                       ;    nprevrec = sf_dyn(jf_tem)%nrec_a(2) 
     126      ENDIF 
     127      ! 
     128      CALL fld_read( kt, 1, sf_dyn )      !=  read data at kt time step   ==! 
     129      ! 
     130      IF( lk_ldfslp .AND. .NOT.lk_c1d )   CALL  dta_dyn_slp( kt )    ! Computation of slopes 
    243131      ! 
    244132      tsn(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:)  * tmask(:,:,:)    ! temperature 
    245133      tsn(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:)  * tmask(:,:,:)    ! salinity 
    246       ! 
     134      wndm(:,:)         = sf_dyn(jf_wnd)%fnow(:,:,1)  * tmask(:,:,1)    ! wind speed - needed for gas exchange 
     135      fmmflx(:,:)       = sf_dyn(jf_fmf)%fnow(:,:,1)  * tmask(:,:,1)    ! downward salt flux (v3.5+) 
     136      fr_i(:,:)         = sf_dyn(jf_ice)%fnow(:,:,1)  * tmask(:,:,1)    ! Sea-ice fraction 
     137      qsr (:,:)         = sf_dyn(jf_qsr)%fnow(:,:,1)  * tmask(:,:,1)    ! solar radiation 
     138      emp (:,:)         = sf_dyn(jf_emp)%fnow(:,:,1)  * tmask(:,:,1)    ! E-P 
     139      IF( ln_dynrnf ) THEN  
     140         rnf (:,:)      = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask(:,:,1)    ! E-P 
     141         IF( ln_dynrnf_depth .AND. lk_vvl )    CALL  dta_dyn_hrnf 
     142      ENDIF 
     143      ! 
     144      un(:,:,:)        = sf_dyn(jf_uwd)%fnow(:,:,:) * umask(:,:,:)    ! effective u-transport 
     145      vn(:,:,:)        = sf_dyn(jf_vwd)%fnow(:,:,:) * vmask(:,:,:)    ! effective v-transport 
     146      wn(:,:,:)        = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:)    ! effective v-transport 
     147      ! 
     148      IF( lk_vvl ) THEN 
     149         CALL wrk_alloc(jpi, jpj, zemp ) 
     150         zhdivtr(:,:,:)    = sf_dyn(jf_div)%fnow(:,:,:)  * tmask(:,:,:)    ! effective u-transport 
     151         emp_b (:,:)       = sf_dyn(jf_empb)%fnow(:,:,1) * tmask(:,:,1)    ! E-P 
     152         zemp(:,:) = 0.5_wp * ( emp(:,:) + emp_b(:,:) ) + rnf(:,:) + fwbcorr * tmask(:,:,1) 
     153         CALL dta_dyn_ssh( kt, zhdivtr, sshb, zemp, ssha, fse3t_a(:,:,:) )  !=  ssh, vertical scale factor & vertical transport 
     154         CALL wrk_dealloc(jpi, jpj, zemp ) 
     155         !                                           Write in the tracer restart file 
     156         !                                          ******************************* 
     157         IF( lrst_trc ) THEN 
     158            IF(lwp) WRITE(numout,*) 
     159            IF(lwp) WRITE(numout,*) 'dta_dyn_ssh : ssh field written in tracer restart file ',   & 
     160               &                    'at it= ', kt,' date= ', ndastp 
     161            IF(lwp) WRITE(numout,*) '~~~~' 
     162            CALL iom_rstput( kt, nitrst, numrtw, 'sshn', ssha ) 
     163            CALL iom_rstput( kt, nitrst, numrtw, 'sshb', sshn ) 
     164         ENDIF 
     165      ENDIF 
    247166      ! 
    248167      CALL eos    ( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 
     
    251170 
    252171      rn2b(:,:,:) = rn2(:,:,:)         ! need for zdfmxl 
    253       CALL zdf_mxl( kt )                                                   ! In any case, we need mxl  
    254       ! 
    255       avt(:,:,:)       = sf_dyn(jf_avt)%fnow(:,:,:)  * tmask(:,:,:)    ! vertical diffusive coefficient  
    256       un (:,:,:)       = sf_dyn(jf_uwd)%fnow(:,:,:)  * umask(:,:,:)    ! u-velocity 
    257       vn (:,:,:)       = sf_dyn(jf_vwd)%fnow(:,:,:)  * vmask(:,:,:)    ! v-velocity  
    258       IF( .NOT.ln_dynwzv ) &                                          ! w-velocity read in file  
    259          wn (:,:,:)    = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:)     
    260       hmld(:,:)        = sf_dyn(jf_mld)%fnow(:,:,1) * tmask(:,:,1)    ! mixed layer depht 
    261       wndm(:,:)        = sf_dyn(jf_wnd)%fnow(:,:,1) * tmask(:,:,1)    ! wind speed - needed for gas exchange 
    262       emp (:,:)        = sf_dyn(jf_emp)%fnow(:,:,1) * tmask(:,:,1)    ! E-P 
    263       fmmflx(:,:)      = sf_dyn(jf_fmf)%fnow(:,:,1) * tmask(:,:,1)    ! downward salt flux (v3.5+) 
    264       fr_i(:,:)        = sf_dyn(jf_ice)%fnow(:,:,1) * tmask(:,:,1)    ! Sea-ice fraction 
    265       qsr (:,:)        = sf_dyn(jf_qsr)%fnow(:,:,1) * tmask(:,:,1)    ! solar radiation 
    266       IF( ln_dynrnf ) & 
    267       rnf (:,:)        = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask(:,:,1)    ! river runoffs  
    268  
    269       !                                                      ! bbl diffusive coef 
     172      CALL zdf_mxl( kt )                                                   ! In any case, we need mxl 
     173      ! 
     174      hmld(:,:)         = sf_dyn(jf_mld)%fnow(:,:,1) * tmask(:,:,1)    ! mixed layer depht 
     175      avt(:,:,:)        = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:)    ! vertical diffusive coefficient  
     176      ! 
    270177#if defined key_trabbl && ! defined key_c1d 
    271       IF( ln_dynbbl ) THEN                                        ! read in a file 
    272          ahu_bbl(:,:)  = sf_dyn(jf_ubl)%fnow(:,:,1) * umask(:,:,1) 
    273          ahv_bbl(:,:)  = sf_dyn(jf_vbl)%fnow(:,:,1) * vmask(:,:,1) 
    274       ELSE                                                        ! Compute bbl coefficients if needed 
    275          tsb(:,:,:,:) = tsn(:,:,:,:) 
    276          CALL bbl( kt, nit000, 'TRC') 
    277       END IF 
     178      ahu_bbl(:,:)      = sf_dyn(jf_ubl)%fnow(:,:,1) * umask(:,:,1)    ! bbl diffusive coef 
     179      ahv_bbl(:,:)      = sf_dyn(jf_vbl)%fnow(:,:,1) * vmask(:,:,1) 
    278180#endif 
    279 #if ( ! defined key_degrad && defined key_traldf_c2d && defined key_traldf_eiv ) && ! defined key_c1d  
    280       aeiw(:,:)        = sf_dyn(jf_eiw)%fnow(:,:,1) * tmask(:,:,1)    ! w-eiv 
    281       !                                                           ! Computes the horizontal values from the vertical value 
    282       DO jj = 2, jpjm1 
    283          DO ji = fs_2, fs_jpim1   ! vector opt. 
    284             aeiu(ji,jj) = .5 * ( aeiw(ji,jj) + aeiw(ji+1,jj  ) )  ! Average the diffusive coefficient at u- v- points 
    285             aeiv(ji,jj) = .5 * ( aeiw(ji,jj) + aeiw(ji  ,jj+1) )  ! at u- v- points 
    286          END DO 
    287       END DO 
    288       CALL lbc_lnk( aeiu, 'U', 1. )   ;   CALL lbc_lnk( aeiv, 'V', 1. )    ! lateral boundary condition 
    289 #endif 
    290        
    291 #if defined key_degrad && ! defined key_c1d  
    292       !                                          ! degrad option : diffusive and eiv coef are 3D 
    293       ahtu(:,:,:) = sf_dyn(jf_ahu)%fnow(:,:,:) * umask(:,:,:) 
    294       ahtv(:,:,:) = sf_dyn(jf_ahv)%fnow(:,:,:) * vmask(:,:,:) 
    295       ahtw(:,:,:) = sf_dyn(jf_ahw)%fnow(:,:,:) * tmask(:,:,:) 
    296 #  if defined key_traldf_eiv  
    297       aeiu(:,:,:) = sf_dyn(jf_eiu)%fnow(:,:,:) * umask(:,:,:) 
    298       aeiv(:,:,:) = sf_dyn(jf_eiv)%fnow(:,:,:) * vmask(:,:,:) 
    299       aeiw(:,:,:) = sf_dyn(jf_eiw)%fnow(:,:,:) * tmask(:,:,:) 
    300 #  endif 
    301 #endif 
     181      ! 
     182      ! 
     183      CALL eos( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 
    302184      ! 
    303185      IF(ln_ctl) THEN                  ! print control 
     
    308190         CALL prt_ctl(tab3d_1=wn               , clinfo1=' wn      - : ', mask1=tmask, ovlap=1, kdim=jpk   ) 
    309191         CALL prt_ctl(tab3d_1=avt              , clinfo1=' kz      - : ', mask1=tmask, ovlap=1, kdim=jpk   ) 
    310          CALL prt_ctl(tab2d_1=fr_i             , clinfo1=' fr_i    - : ', mask1=tmask, ovlap=1 ) 
    311          CALL prt_ctl(tab2d_1=hmld             , clinfo1=' hmld    - : ', mask1=tmask, ovlap=1 ) 
    312          CALL prt_ctl(tab2d_1=fmmflx           , clinfo1=' fmmflx  - : ', mask1=tmask, ovlap=1 ) 
    313          CALL prt_ctl(tab2d_1=emp              , clinfo1=' emp     - : ', mask1=tmask, ovlap=1 ) 
    314          CALL prt_ctl(tab2d_1=wndm             , clinfo1=' wspd    - : ', mask1=tmask, ovlap=1 ) 
    315          CALL prt_ctl(tab2d_1=qsr              , clinfo1=' qsr     - : ', mask1=tmask, ovlap=1 ) 
     192!         CALL prt_ctl(tab2d_1=fr_i             , clinfo1=' fr_i    - : ', mask1=tmask, ovlap=1 ) 
     193!         CALL prt_ctl(tab2d_1=hmld             , clinfo1=' hmld    - : ', mask1=tmask, ovlap=1 ) 
     194!         CALL prt_ctl(tab2d_1=fmmflx           , clinfo1=' fmmflx  - : ', mask1=tmask, ovlap=1 ) 
     195!         CALL prt_ctl(tab2d_1=emp              , clinfo1=' emp     - : ', mask1=tmask, ovlap=1 ) 
     196!         CALL prt_ctl(tab2d_1=wndm             , clinfo1=' wspd    - : ', mask1=tmask, ovlap=1 ) 
     197!         CALL prt_ctl(tab2d_1=qsr              , clinfo1=' qsr     - : ', mask1=tmask, ovlap=1 ) 
    316198      ENDIF 
    317199      ! 
     
    335217      INTEGER  :: inum, idv, idimv                   ! local integer 
    336218      INTEGER  :: ios                                ! Local integer output status for namelist read 
    337       !! 
    338       CHARACTER(len=100)            ::  cn_dir   !   Root directory for location of core files 
    339       TYPE(FLD_N), DIMENSION(jpfld) ::  slf_d    ! array of namelist informations on the fields to read 
    340       TYPE(FLD_N) :: sn_tem, sn_sal, sn_mld, sn_emp, sn_ice, sn_qsr, sn_wnd, sn_rnf  ! informations about the fields to be read 
    341       TYPE(FLD_N) :: sn_uwd, sn_vwd, sn_wwd, sn_avt, sn_ubl, sn_vbl          !   "                                 " 
    342       TYPE(FLD_N) :: sn_ahu, sn_ahv, sn_ahw, sn_eiu, sn_eiv, sn_eiw, sn_fmf  !   "                                 " 
    343       !!---------------------------------------------------------------------- 
    344       ! 
    345       NAMELIST/namdta_dyn/cn_dir, ln_dynwzv, ln_dynbbl, ln_degrad, ln_dynrnf,    & 
    346          &                sn_tem, sn_sal, sn_mld, sn_emp, sn_ice, sn_qsr, sn_wnd, sn_rnf,  & 
    347          &                sn_uwd, sn_vwd, sn_wwd, sn_avt, sn_ubl, sn_vbl,          & 
    348          &                sn_ahu, sn_ahv, sn_ahw, sn_eiu, sn_eiv, sn_eiw, sn_fmf 
     219      INTEGER  :: ji, jj, jk 
     220      REAL(wp) :: zcoef 
     221      INTEGER  :: nkrnf_max 
     222      REAL(wp) :: hrnf_max 
     223      !! 
     224      CHARACTER(len=100)            ::  cn_dir        !   Root directory for location of core files 
     225      TYPE(FLD_N), DIMENSION(jpfld) ::  slf_d         ! array of namelist informations on the fields to read 
     226      TYPE(FLD_N) :: sn_uwd, sn_vwd, sn_wwd, sn_empb, sn_emp  ! informations about the fields to be read 
     227      TYPE(FLD_N) :: sn_tem , sn_sal , sn_avt   !   "                 " 
     228      TYPE(FLD_N) :: sn_mld, sn_qsr, sn_wnd , sn_ice , sn_fmf   !   "               " 
     229      TYPE(FLD_N) :: sn_ubl, sn_vbl, sn_rnf    !   "              " 
     230      TYPE(FLD_N) :: sn_div  ! informations about the fields to be read 
     231      !!---------------------------------------------------------------------- 
     232 
     233      NAMELIST/namdta_dyn/cn_dir, ln_dynrnf, ln_dynrnf_depth,  fwbcorr, & 
     234         &                sn_uwd, sn_vwd, sn_wwd, sn_emp,    & 
     235         &                sn_avt, sn_tem, sn_sal, sn_mld , sn_qsr ,   & 
     236         &                sn_wnd, sn_ice, sn_fmf,                    & 
     237         &                sn_ubl, sn_vbl, sn_rnf,                   & 
     238         &                sn_empb, sn_div  
    349239      ! 
    350240      REWIND( numnam_ref )              ! Namelist namdta_dyn in reference namelist : Offline: init. of dynamical data 
     
    363253         WRITE(numout,*) '~~~~~~~ ' 
    364254         WRITE(numout,*) '   Namelist namdta_dyn' 
    365          WRITE(numout,*) '      vertical velocity read from file (T) or computed (F) ln_dynwzv  = ', ln_dynwzv 
    366          WRITE(numout,*) '      bbl coef read from file (T) or computed (F)          ln_dynbbl  = ', ln_dynbbl 
    367          WRITE(numout,*) '      degradation option enabled (T) or not (F)            ln_degrad  = ', ln_degrad 
    368          WRITE(numout,*) '      river runoff option enabled (T) or not (F)           ln_dynrnf  = ', ln_dynrnf 
     255         WRITE(numout,*) '      runoffs option enabled (T) or not (F)            ln_dynrnf        = ', ln_dynrnf 
     256         WRITE(numout,*) '      runoffs is spread in vertical                    ln_dynrnf_depth  = ', ln_dynrnf_depth 
     257         WRITE(numout,*) '      annual global mean of empmr for ssh correction   fwbcorr          = ', fwbcorr 
    369258         WRITE(numout,*) 
    370259      ENDIF 
    371260      !  
    372       IF( ln_degrad .AND. .NOT.lk_degrad ) THEN 
    373          CALL ctl_warn( 'dta_dyn_init: degradation option requires key_degrad activated ; force ln_degrad to false' ) 
    374          ln_degrad = .FALSE. 
    375       ENDIF 
    376       IF( ln_dynbbl .AND. ( .NOT.lk_trabbl .OR. lk_c1d ) ) THEN 
    377          CALL ctl_warn( 'dta_dyn_init: bbl option requires key_trabbl activated ; force ln_dynbbl to false' ) 
    378          ln_dynbbl = .FALSE. 
    379       ENDIF 
    380  
    381       jf_tem = 1   ;   jf_sal = 2   ;  jf_mld = 3   ;  jf_emp = 4   ;   jf_fmf  = 5   ;  jf_ice = 6   ;   jf_qsr = 7 
    382       jf_wnd = 8   ;   jf_uwd = 9   ;  jf_vwd = 10  ;  jf_wwd = 11  ;   jf_avt  = 12  ;  jfld  = jf_avt 
    383       ! 
    384       slf_d(jf_tem) = sn_tem   ;   slf_d(jf_sal)  = sn_sal   ;   slf_d(jf_mld) = sn_mld 
    385       slf_d(jf_emp) = sn_emp   ;   slf_d(jf_fmf ) = sn_fmf   ;   slf_d(jf_ice) = sn_ice  
    386       slf_d(jf_qsr) = sn_qsr   ;   slf_d(jf_wnd)  = sn_wnd   ;   slf_d(jf_avt) = sn_avt  
    387       slf_d(jf_uwd) = sn_uwd   ;   slf_d(jf_vwd)  = sn_vwd   ;   slf_d(jf_wwd) = sn_wwd 
    388  
     261 
     262      jf_uwd  = 1     ;   jf_vwd  = 2    ;   jf_wwd = 3    ;   jf_emp = 4    ;   jf_avt = 5 
     263      jf_tem  = 6     ;   jf_sal  = 7    ;   jf_mld = 8    ;   jf_qsr = 9 
     264      jf_wnd  = 10    ;   jf_ice  = 11   ;   jf_fmf = 12   ;   jfld   = jf_fmf 
     265 
     266      ! 
     267      slf_d(jf_uwd)  = sn_uwd    ;   slf_d(jf_vwd)  = sn_vwd   ;   slf_d(jf_wwd) = sn_wwd 
     268      slf_d(jf_emp)  = sn_emp    ;   slf_d(jf_avt)  = sn_avt 
     269      slf_d(jf_tem)  = sn_tem    ;   slf_d(jf_sal)  = sn_sal   ;   slf_d(jf_mld) = sn_mld 
     270      slf_d(jf_qsr)  = sn_qsr    ;   slf_d(jf_wnd)  = sn_wnd   ;   slf_d(jf_ice) = sn_ice 
     271      slf_d(jf_fmf)  = sn_fmf 
     272 
     273 
     274      ! 
     275      IF( lk_vvl ) THEN 
     276                 jf_div  = jfld + 1    ;         jf_empb  = jfld + 2      ;      jfld = jf_empb 
     277           slf_d(jf_div) = sn_div      ;   slf_d(jf_empb) = sn_empb 
     278      ENDIF 
     279      ! 
     280      IF( lk_trabbl ) THEN 
     281                 jf_ubl  = jfld + 1    ;         jf_vbl  = jfld + 2     ;      jfld = jf_vbl 
     282           slf_d(jf_ubl) = sn_ubl      ;   slf_d(jf_vbl) = sn_vbl 
     283      ENDIF 
    389284      ! 
    390285      IF( ln_dynrnf ) THEN 
    391                 jf_rnf = jfld + 1  ;  jfld  = jf_rnf 
    392          slf_d(jf_rnf) = sn_rnf 
     286                jf_rnf  = jfld + 1     ;     jfld  = jf_rnf 
     287          slf_d(jf_rnf) = sn_rnf     
    393288      ELSE 
    394          rnf (:,:) = 0._wp 
    395       ENDIF 
    396  
    397       ! 
    398       IF( .NOT.ln_degrad ) THEN     ! no degrad option 
    399          IF( lk_traldf_eiv .AND. ln_dynbbl ) THEN        ! eiv & bbl 
    400                  jf_ubl  = jfld + 1 ;        jf_vbl  = jfld + 2 ;        jf_eiw  = jfld + 3   ;   jfld = jf_eiw 
    401            slf_d(jf_ubl) = sn_ubl  ;   slf_d(jf_vbl) = sn_vbl  ;   slf_d(jf_eiw) = sn_eiw 
    402          ENDIF 
    403          IF( .NOT.lk_traldf_eiv .AND. ln_dynbbl ) THEN   ! no eiv & bbl 
    404                  jf_ubl  = jfld + 1 ;        jf_vbl  = jfld + 2 ;  jfld = jf_vbl 
    405            slf_d(jf_ubl) = sn_ubl  ;   slf_d(jf_vbl) = sn_vbl 
    406          ENDIF 
    407          IF( lk_traldf_eiv .AND. .NOT.ln_dynbbl ) THEN   ! eiv & no bbl 
    408            jf_eiw = jfld + 1 ; jfld = jf_eiw ; slf_d(jf_eiw) = sn_eiw 
    409          ENDIF 
    410       ELSE 
    411               jf_ahu  = jfld + 1 ;        jf_ahv  = jfld + 2 ;        jf_ahw  = jfld + 3  ;  jfld = jf_ahw 
    412         slf_d(jf_ahu) = sn_ahu  ;   slf_d(jf_ahv) = sn_ahv  ;   slf_d(jf_ahw) = sn_ahw 
    413         IF( lk_traldf_eiv .AND. ln_dynbbl ) THEN         ! eiv & bbl 
    414                  jf_ubl  = jfld + 1 ;        jf_vbl  = jfld + 2 ; 
    415            slf_d(jf_ubl) = sn_ubl  ;   slf_d(jf_vbl) = sn_vbl 
    416                  jf_eiu  = jfld + 3 ;        jf_eiv  = jfld + 4 ;    jf_eiw  = jfld + 5   ;  jfld = jf_eiw  
    417            slf_d(jf_eiu) = sn_eiu  ;   slf_d(jf_eiv) = sn_eiv  ;    slf_d(jf_eiw) = sn_eiw 
    418         ENDIF 
    419         IF( .NOT.lk_traldf_eiv .AND. ln_dynbbl ) THEN    ! no eiv & bbl 
    420                  jf_ubl  = jfld + 1 ;        jf_vbl  = jfld + 2 ;  jfld = jf_vbl 
    421            slf_d(jf_ubl) = sn_ubl  ;   slf_d(jf_vbl) = sn_vbl 
    422         ENDIF 
    423         IF( lk_traldf_eiv .AND. .NOT.ln_dynbbl ) THEN    ! eiv & no bbl 
    424                  jf_eiu  = jfld + 1 ;         jf_eiv  = jfld + 2 ;    jf_eiw  = jfld + 3   ; jfld = jf_eiw  
    425            slf_d(jf_eiu) = sn_eiu  ;   slf_d(jf_eiv) = sn_eiv  ;   slf_d(jf_eiw) = sn_eiw 
    426         ENDIF 
    427       ENDIF 
     289         rnf(:,:) = 0._wp 
     290      ENDIF 
     291 
    428292   
    429293      ALLOCATE( sf_dyn(jfld), STAT=ierr )         ! set sf structure 
    430       IF( ierr > 0 ) THEN 
     294      IF( ierr > 0 )  THEN 
    431295         CALL ctl_stop( 'dta_dyn: unable to allocate sf structure' )   ;   RETURN 
    432296      ENDIF 
    433297      !                                         ! fill sf with slf_i and control print 
    434298      CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' ) 
     299      ! 
    435300      ! Open file for each variable to get his number of dimension 
    436301      DO ifpr = 1, jfld 
     
    456321            ALLOCATE( uslpdta (jpi,jpj,jpk,2), vslpdta (jpi,jpj,jpk,2),    & 
    457322            &         wslpidta(jpi,jpj,jpk,2), wslpjdta(jpi,jpj,jpk,2), STAT=ierr2 ) 
    458          ELSE 
    459             ALLOCATE( uslpnow (jpi,jpj,jpk)  , vslpnow (jpi,jpj,jpk)  ,    & 
    460             &         wslpinow(jpi,jpj,jpk)  , wslpjnow(jpi,jpj,jpk)  , STAT=ierr2 ) 
    461          ENDIF  
    462          IF( ierr2 > 0 ) THEN 
    463             CALL ctl_stop( 'dta_dyn_init : unable to allocate slope arrays' )   ;   RETURN 
     323            ! 
     324            IF( ierr2 > 0 )  THEN 
     325               CALL ctl_stop( 'dta_dyn_init : unable to allocate slope arrays' )   ;   RETURN 
     326            ENDIF 
    464327         ENDIF 
    465328      ENDIF 
    466       IF( ln_dynwzv ) THEN                  ! slopes  
    467          IF( sf_dyn(jf_uwd)%ln_tint ) THEN      ! time interpolation 
    468             ALLOCATE( wdta(jpi,jpj,jpk,2), STAT=ierr3 ) 
    469          ELSE 
    470             ALLOCATE( wnow(jpi,jpj,jpk)  , STAT=ierr3 ) 
    471          ENDIF  
    472          IF( ierr3 > 0 ) THEN 
    473             CALL ctl_stop( 'dta_dyn_init : unable to allocate wdta arrays' )   ;   RETURN 
    474          ENDIF 
    475       ENDIF 
    476       ! 
    477       CALL dta_dyn( nit000 ) 
    478       ! 
    479    END SUBROUTINE dta_dyn_init 
    480  
    481    SUBROUTINE dta_dyn_wzv( pu, pv, pw ) 
    482       !!---------------------------------------------------------------------- 
    483       !!                    ***  ROUTINE wzv  *** 
    484       !! 
    485       !! ** Purpose :   Compute the now vertical velocity after the array swap 
    486       !! 
    487       !! ** Method  : - compute the now divergence given by : 
    488       !!         * z-coordinate ONLY !!!! 
    489       !!         hdiv = 1/(e1t*e2t) [ di(e2u  u) + dj(e1v  v) ] 
    490       !!     - Using the incompressibility hypothesis, the vertical 
    491       !!      velocity is computed by integrating the horizontal divergence 
    492       !!      from the bottom to the surface. 
    493       !!        The boundary conditions are w=0 at the bottom (no flux). 
    494       !!---------------------------------------------------------------------- 
    495       USE oce, ONLY:  zhdiv => hdivn 
    496       ! 
    497       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pu, pv    !:  horizontal velocities 
    498       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) :: pw        !:  vertical velocity 
    499       !! 
    500       INTEGER  ::  ji, jj, jk 
    501       REAL(wp) ::  zu, zu1, zv, zv1, zet 
    502       !!---------------------------------------------------------------------- 
    503       ! 
    504       ! Computation of vertical velocity using horizontal divergence 
    505       zhdiv(:,:,:) = 0._wp 
    506       DO jk = 1, jpkm1 
    507          DO jj = 2, jpjm1 
    508             DO ji = fs_2, fs_jpim1   ! vector opt. 
    509                zu  = pu(ji  ,jj  ,jk) * umask(ji  ,jj  ,jk) * e2u(ji  ,jj  ) * fse3u(ji  ,jj  ,jk) 
    510                zu1 = pu(ji-1,jj  ,jk) * umask(ji-1,jj  ,jk) * e2u(ji-1,jj  ) * fse3u(ji-1,jj  ,jk) 
    511                zv  = pv(ji  ,jj  ,jk) * vmask(ji  ,jj  ,jk) * e1v(ji  ,jj  ) * fse3v(ji  ,jj  ,jk) 
    512                zv1 = pv(ji  ,jj-1,jk) * vmask(ji  ,jj-1,jk) * e1v(ji  ,jj-1) * fse3v(ji  ,jj-1,jk) 
    513                zet = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    514                zhdiv(ji,jj,jk) = ( zu - zu1 + zv - zv1 ) * zet  
     329      ! 
     330      IF( lk_vvl ) THEN 
     331        IF( .NOT. sf_dyn(jf_uwd)%ln_clim .AND. ln_rsttr .AND.    &                     ! Restart: read in restart file 
     332           iom_varid( numrtr, 'sshn', ldstop = .FALSE. ) > 0 ) THEN 
     333           IF(lwp) WRITE(numout,*) ' sshn forcing fields read in the restart file for initialisation' 
     334           CALL iom_get( numrtr, jpdom_autoglo, 'sshn', sshn(:,:)   ) 
     335           CALL iom_get( numrtr, jpdom_autoglo, 'sshb', sshb(:,:)   ) 
     336        ELSE 
     337           IF(lwp) WRITE(numout,*) ' sshn forcing fields read in the restart file for initialisation' 
     338           CALL iom_open( 'restart', inum ) 
     339           CALL iom_get( inum, jpdom_autoglo, 'sshn', sshn(:,:)   ) 
     340           CALL iom_get( inum, jpdom_autoglo, 'sshb', sshb(:,:)   ) 
     341           CALL iom_close( inum )                                        ! close file 
     342        ENDIF 
     343        ! 
     344        DO jk = 1, jpkm1 
     345           fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) ) 
     346        ENDDO 
     347        fse3t_a(:,:,jpk) = e3t_0(:,:,jpk) 
     348 
     349        ! Horizontal scale factor interpolations 
     350        ! -------------------------------------- 
     351        CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' ) 
     352        CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 
     353 
     354        ! Vertical scale factor interpolations 
     355        ! ------------------------------------ 
     356        CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n(:,:,:), 'W' ) 
     357   
     358        fse3t_b(:,:,:)  = fse3t_n(:,:,:) 
     359        fse3u_b(:,:,:)  = fse3u_n(:,:,:) 
     360        fse3v_b(:,:,:)  = fse3v_n(:,:,:) 
     361 
     362        ! t- and w- points depth 
     363        ! ---------------------- 
     364        fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 
     365        fsdepw_n(:,:,1) = 0.0_wp 
     366 
     367        DO jk = 2, jpk 
     368           DO jj = 1,jpj 
     369              DO ji = 1,jpi 
     370                !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere 
     371                !    tmask = wmask, ie everywhere expect at jk = mikt 
     372                                                                   ! 1 for jk = 
     373                                                                   ! mikt 
     374                 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
     375                 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 
     376                 fsdept_n(ji,jj,jk) =      zcoef  * ( fsdepw_n(ji,jj,jk  ) + 0.5 * fse3w_n(ji,jj,jk))  & 
     377                     &                + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)) 
     378              END DO 
     379           END DO 
     380        END DO 
     381 
     382        fsdept_b(:,:,:) = fsdept_n(:,:,:) 
     383        fsdepw_b(:,:,:) = fsdepw_n(:,:,:) 
     384        ! 
     385      ENDIF 
     386      ! 
     387      IF( ln_dynrnf .AND. ln_dynrnf_depth ) THEN       ! read depht over which runoffs are distributed 
     388         IF(lwp) WRITE(numout,*)  
     389         IF(lwp) WRITE(numout,*) ' read in the file depht over which runoffs are distributed' 
     390         CALL iom_open ( "runoffs", inum )                           ! open file 
     391         CALL iom_get  ( inum, jpdom_data, 'rodepth', h_rnf )   ! read the river mouth array 
     392         CALL iom_close( inum )                                        ! close file 
     393         ! 
     394         nk_rnf(:,:) = 0                               ! set the number of level over which river runoffs are applied 
     395         DO jj = 1, jpj 
     396            DO ji = 1, jpi 
     397               IF( h_rnf(ji,jj) > 0._wp ) THEN 
     398                  jk = 2 
     399                  DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1 
     400                  END DO 
     401                  nk_rnf(ji,jj) = jk 
     402               ELSEIF( h_rnf(ji,jj) == -1._wp   ) THEN   ;  nk_rnf(ji,jj) = 1 
     403               ELSEIF( h_rnf(ji,jj) == -999._wp ) THEN   ;  nk_rnf(ji,jj) = mbkt(ji,jj) 
     404               ELSE 
     405                  CALL ctl_stop( 'sbc_rnf_init: runoff depth not positive, and not -999 or -1, rnf value in file fort.999'  ) 
     406                  WRITE(999,*) 'ji, jj, h_rnf(ji,jj) :', ji, jj, h_rnf(ji,jj) 
     407               ENDIF 
    515408            END DO 
    516409         END DO 
     410         DO jj = 1, jpj                                ! set the associated depth 
     411            DO ji = 1, jpi 
     412               h_rnf(ji,jj) = 0._wp 
     413               DO jk = 1, nk_rnf(ji,jj) 
     414                  h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk) 
     415               END DO 
     416            END DO 
     417         END DO 
     418      ELSE                                       ! runoffs applied at the surface 
     419         nk_rnf(:,:) = 1 
     420         h_rnf (:,:) = fse3t(:,:,1) 
     421      ENDIF 
     422      nkrnf_max = MAXVAL( nk_rnf(:,:) ) 
     423      hrnf_max = MAXVAL( h_rnf(:,:) ) 
     424      IF( lk_mpp )  THEN 
     425         CALL mpp_max( nkrnf_max )                 ! max over the  global domain 
     426         CALL mpp_max( hrnf_max )                 ! max over the  global domain 
     427      ENDIF 
     428      IF(lwp) WRITE(numout,*) ' ' 
     429      IF(lwp) WRITE(numout,*) ' max depht of runoff : ', hrnf_max,'    max level  : ', nkrnf_max 
     430      IF(lwp) WRITE(numout,*) ' ' 
     431      ! 
     432      CALL dta_dyn( nit000 ) 
     433      ! 
     434   END SUBROUTINE dta_dyn_init 
     435 
     436   SUBROUTINE dta_dyn_swp( kt ) 
     437     !!--------------------------------------------------------------------- 
     438      !!                    ***  ROUTINE dta_dyn_swp  *** 
     439      !! 
     440      !! ** Purpose : Swap and the data and compute the vertical scale factor at U/V/W point 
     441      !!              and the depht 
     442      !! 
     443      !!--------------------------------------------------------------------- 
     444      INTEGER, INTENT(in) :: kt       ! time step 
     445      INTEGER             :: ji, jj, jk 
     446      REAL(wp)            :: zcoef 
     447      ! 
     448      !!--------------------------------------------------------------------- 
     449 
     450      IF( kt == nit000 ) THEN 
     451         IF(lwp) WRITE(numout,*) 
     452         IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height' 
     453         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
     454      ENDIF 
     455 
     456      sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:))  ! before <-- now filtered 
     457      sshn(:,:) = ssha(:,:) 
     458 
     459      fse3t_n(:,:,:) = fse3t_a(:,:,:) 
     460 
     461      ! Reconstruction of all vertical scale factors at now and before time steps 
     462      ! ============================================================================= 
     463 
     464      ! Horizontal scale factor interpolations 
     465      ! -------------------------------------- 
     466      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' ) 
     467      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 
     468 
     469      ! Vertical scale factor interpolations 
     470      ! ------------------------------------ 
     471      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' ) 
     472 
     473      fse3t_b(:,:,:)  = fse3t_n(:,:,:) 
     474      fse3u_b(:,:,:)  = fse3u_n(:,:,:) 
     475      fse3v_b(:,:,:)  = fse3v_n(:,:,:) 
     476 
     477      ! t- and w- points depth 
     478      ! ---------------------- 
     479      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 
     480      fsdepw_n(:,:,1) = 0.0_wp 
     481 
     482      DO jk = 2, jpk 
     483         DO jj = 1,jpj 
     484            DO ji = 1,jpi 
     485                 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
     486                 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 
     487                 fsdept_n(ji,jj,jk) =      zcoef  * ( fsdepw_n(ji,jj,jk  ) + 0.5 * fse3w_n(ji,jj,jk))  & 
     488                     &                + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)) 
     489              END DO 
     490           END DO 
     491        END DO 
     492 
     493      fsdept_b(:,:,:) = fsdept_n(:,:,:) 
     494      fsdepw_b(:,:,:) = fsdepw_n(:,:,:) 
     495 
     496      ! 
     497   END SUBROUTINE dta_dyn_swp 
     498 
     499   SUBROUTINE dta_dyn_ssh( kt, phdivtr, psshb,  pemp, pssha, pe3ta ) 
     500      !!---------------------------------------------------------------------- 
     501      !!                ***  ROUTINE dta_dyn_wzv  *** 
     502      !!                    
     503      !! ** Purpose :   compute the after ssh (ssha) and the now vertical velocity 
     504      !! 
     505      !! ** Method  : Using the incompressibility hypothesis,  
     506      !!        - the ssh increment is computed by integrating the horizontal divergence  
     507      !!          and multiply by the time step. 
     508      !! 
     509      !!        - compute the after scale factor : repartition of ssh INCREMENT proportionnaly 
     510      !!                                           to the level thickness ( z-star case ) 
     511      !! 
     512      !!        - the vertical velocity is computed by integrating the horizontal divergence   
     513      !!          from the bottom to the surface minus the scale factor evolution. 
     514      !!          The boundary conditions are w=0 at the bottom (no flux) 
     515      !! 
     516      !! ** action  :   ssha / e3t_a / wn 
     517      !! 
     518      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
     519      !!---------------------------------------------------------------------- 
     520      !! * Arguments 
     521      INTEGER,                                   INTENT(in )    :: kt        !  time-step 
     522      REAL(wp), DIMENSION(jpi,jpj,jpk)          , INTENT(in )   :: phdivtr   ! horizontal divergence transport 
     523      REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(in )   :: psshb     ! now ssh 
     524      REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(in )   :: pemp      ! evaporation minus precipitation 
     525      REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(inout) :: pssha     ! after ssh 
     526      REAL(wp), DIMENSION(jpi,jpj,jpk), OPTIONAL, INTENT(out)   :: pe3ta     ! after vertical scale factor 
     527      !! * Local declarations 
     528      INTEGER                       :: jk 
     529      REAL(wp), DIMENSION(jpi,jpj)  :: zhdiv   
     530      REAL(wp)  :: z2dt   
     531      !!---------------------------------------------------------------------- 
     532       
     533      ! 
     534      z2dt = 2._wp * rdt 
     535      ! 
     536      zhdiv(:,:) = 0._wp 
     537      DO jk = 1, jpkm1 
     538         zhdiv(:,:) = zhdiv(:,:) +  phdivtr(:,:,jk) * tmask(:,:,jk) 
    517539      END DO 
    518       CALL lbc_lnk( zhdiv, 'T', 1. )      ! Lateral boundary conditions on zhdiv 
    519       ! 
    520       ! computation of vertical velocity from the bottom 
    521       pw(:,:,jpk) = 0._wp 
    522       DO jk = jpkm1, 1, -1 
    523          pw(:,:,jk) = pw(:,:,jk+1) - fse3t(:,:,jk) * zhdiv(:,:,jk) 
     540      !                                                ! Sea surface  elevation time-stepping 
     541      pssha(:,:) = ( psshb(:,:) - z2dt * ( r1_rau0 * pemp(:,:)  + zhdiv(:,:) ) ) * ssmask(:,:) 
     542      !                                                 !  
     543      !                                                 ! After acale factors at t-points ( z_star coordinate ) 
     544      DO jk = 1, jpkm1 
     545        pe3ta(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + pssha(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) ) 
    524546      END DO 
    525547      ! 
    526    END SUBROUTINE dta_dyn_wzv 
    527  
    528    SUBROUTINE dta_dyn_slp( kt, pts, puslp, pvslp, pwslpi, pwslpj ) 
     548   END SUBROUTINE dta_dyn_ssh 
     549 
     550 
     551   SUBROUTINE dta_dyn_hrnf 
     552      !!---------------------------------------------------------------------- 
     553      !!                  ***  ROUTINE sbc_rnf  *** 
     554      !! 
     555      !! ** Purpose :   update the horizontal divergence with the runoff inflow 
     556      !! 
     557      !! ** Method  : 
     558      !!                CAUTION : rnf is positive (inflow) decreasing the 
     559      !!                          divergence and expressed in m/s 
     560      !! 
     561      !! ** Action  :   phdivn   decreased by the runoff inflow 
     562      !!---------------------------------------------------------------------- 
     563      !! 
     564      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     565      !!---------------------------------------------------------------------- 
     566      ! 
     567      DO jj = 1, jpj                   ! update the depth over which runoffs are distributed 
     568         DO ji = 1, jpi 
     569            h_rnf(ji,jj) = 0._wp 
     570            DO jk = 1, nk_rnf(ji,jj)                           ! recalculates h_rnf to be the depth in metres 
     571                h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk)   ! to the bottom of the relevant grid box 
     572            END DO 
     573        END DO 
     574      END DO 
     575      ! 
     576   END SUBROUTINE dta_dyn_hrnf 
     577 
     578 
     579 
     580   SUBROUTINE dta_dyn_slp( kt ) 
     581      !!--------------------------------------------------------------------- 
     582      !!                    ***  ROUTINE dta_dyn_slp  *** 
     583      !! 
     584      !! ** Purpose : Computation of slope 
     585      !! 
     586      !!--------------------------------------------------------------------- 
     587      USE oce, ONLY:  zts => tsa  
     588      ! 
     589      INTEGER,  INTENT(in) :: kt       ! time step 
     590      ! 
     591      INTEGER  ::   ji, jj     ! dummy loop indices 
     592      REAL(wp) ::   ztinta     ! ratio applied to after  records when doing time interpolation 
     593      REAL(wp) ::   ztintb     ! ratio applied to before records when doing time interpolation 
     594      INTEGER  ::   iswap 
     595      REAL(wp), POINTER, DIMENSION(:,:,:) :: zuslp, zvslp, zwslpi, zwslpj 
     596      !!--------------------------------------------------------------------- 
     597      ! 
     598      CALL wrk_alloc(jpi, jpj, jpk, zuslp, zvslp, zwslpi, zwslpj ) 
     599      ! 
     600      IF( sf_dyn(jf_tem)%ln_tint ) THEN    ! Computes slopes (here avt is used as workspace)                        
     601         IF( kt == nit000 ) THEN 
     602            zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,1) * tmask(:,:,:)   ! temperature 
     603            zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:)   ! salinity  
     604            avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,1) * tmask(:,:,:)   ! vertical diffusive coef. 
     605            CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
     606            uslpdta (:,:,:,1) = zuslp (:,:,:)  
     607            vslpdta (:,:,:,1) = zvslp (:,:,:)  
     608            wslpidta(:,:,:,1) = zwslpi(:,:,:)  
     609            wslpjdta(:,:,:,1) = zwslpj(:,:,:)  
     610            ! 
     611            zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:)   ! temperature 
     612            zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity  
     613            avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef. 
     614            CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
     615            uslpdta (:,:,:,2) = zuslp (:,:,:)  
     616            vslpdta (:,:,:,2) = zvslp (:,:,:)  
     617            wslpidta(:,:,:,2) = zwslpi(:,:,:)  
     618            wslpjdta(:,:,:,2) = zwslpj(:,:,:)  
     619         ELSE 
     620           !  
     621           iswap = 0 
     622           IF( sf_dyn(jf_tem)%nrec_a(2) - nprevrec /= 0 )  iswap = 1 
     623           IF( nsecdyn > sf_dyn(jf_tem)%nrec_b(2) .AND. iswap == 1 )  THEN    ! read/update the after data 
     624              IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt 
     625              uslpdta (:,:,:,1) =  uslpdta (:,:,:,2)         ! swap the data 
     626              vslpdta (:,:,:,1) =  vslpdta (:,:,:,2)   
     627              wslpidta(:,:,:,1) =  wslpidta(:,:,:,2)  
     628              wslpjdta(:,:,:,1) =  wslpjdta(:,:,:,2)  
     629              ! 
     630              zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:)   ! temperature 
     631              zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity  
     632              avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef. 
     633              CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
     634              ! 
     635              uslpdta (:,:,:,2) = zuslp (:,:,:)  
     636              vslpdta (:,:,:,2) = zvslp (:,:,:)  
     637              wslpidta(:,:,:,2) = zwslpi(:,:,:)  
     638              wslpjdta(:,:,:,2) = zwslpj(:,:,:)  
     639            ENDIF 
     640         ENDIF 
     641      ENDIF 
     642      ! 
     643      IF( sf_dyn(jf_tem)%ln_tint )  THEN 
     644         ztinta =  REAL( nsecdyn - sf_dyn(jf_tem)%nrec_b(2), wp )  & 
     645            &    / REAL( sf_dyn(jf_tem)%nrec_a(2) - sf_dyn(jf_tem)%nrec_b(2), wp ) 
     646         ztintb =  1. - ztinta 
     647#if defined key_ldfslp && ! defined key_c1d 
     648         uslp (:,:,:) = ztintb * uslpdta (:,:,:,1)  + ztinta * uslpdta (:,:,:,2)   
     649         vslp (:,:,:) = ztintb * vslpdta (:,:,:,1)  + ztinta * vslpdta (:,:,:,2)   
     650         wslpi(:,:,:) = ztintb * wslpidta(:,:,:,1)  + ztinta * wslpidta(:,:,:,2)   
     651         wslpj(:,:,:) = ztintb * wslpjdta(:,:,:,1)  + ztinta * wslpjdta(:,:,:,2)   
     652#endif 
     653      ELSE 
     654         zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:)   ! temperature 
     655         zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:)   ! salinity  
     656         avt(:,:,:)        = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:)   ! vertical diffusive coef. 
     657         CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
     658         ! 
     659#if defined key_ldfslp && ! defined key_c1d 
     660         uslp (:,:,:) = zuslp (:,:,:) 
     661         vslp (:,:,:) = zvslp (:,:,:) 
     662         wslpi(:,:,:) = zwslpi(:,:,:) 
     663         wslpj(:,:,:) = zwslpj(:,:,:) 
     664#endif 
     665      ENDIF 
     666      ! 
     667      CALL wrk_dealloc(jpi, jpj, jpk, zuslp, zvslp, zwslpi, zwslpj ) 
     668      ! 
     669   END SUBROUTINE dta_dyn_slp 
     670 
     671   SUBROUTINE compute_slopes( kt, pts, puslp, pvslp, pwslpi, pwslpj ) 
    529672      !!--------------------------------------------------------------------- 
    530673      !!                    ***  ROUTINE dta_dyn_slp  *** 
     
    568711#endif 
    569712      ! 
    570    END SUBROUTINE dta_dyn_slp 
     713   END SUBROUTINE compute_slopes 
     714 
    571715   !!====================================================================== 
    572716END MODULE dtadyn 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OFF_SRC/nemogcm.F90

    r5504 r7522  
    3434   USE trcstp          ! passive tracer time-stepping      (trc_stp routine) 
    3535   USE dtadyn          ! Lecture and interpolation of the dynamical fields 
    36    !              ! I/O & MPP 
     36   !                   ! I/O & MPP 
    3737   USE iom             ! I/O library 
    3838   USE in_out_manager  ! I/O manager 
     
    5050   USE trcnam 
    5151   USE trcrst 
    52    USE diaptr         ! Need to initialise this as some variables are used in if statements later 
    5352 
    5453   IMPLICIT NONE 
     
    9493      istp = nit000 
    9594      !  
    96       CALL iom_init( cxios_context )            ! iom_put initialization (must be done after nemo_init for AGRIF+XIOS+OASIS) 
    97       !  
    9895      DO WHILE ( istp <= nitend .AND. nstop == 0 )    ! time stepping 
    9996         ! 
    100          IF( istp /= nit000 )   CALL day      ( istp )         ! Calendar (day was already called at nit000 in day_init) 
    101                                 CALL iom_setkt( istp - nit000 + 1, "nemo" )   ! say to iom that we are at time step kstp 
    102                                 CALL dta_dyn  ( istp )         ! Interpolation of the dynamical fields 
    103                                 CALL trc_stp  ( istp )         ! time-stepping 
    104                                 CALL stp_ctl  ( istp, indic )  ! Time loop: control and print 
     97         IF( istp == nit000 )  CALL iom_init( cxios_context )            ! iom_put initialization 
     98         IF( istp /= nit000 )   CALL day        ( istp )         ! Calendar (day was already called at nit000 in day_init) 
     99                                CALL iom_setkt  ( istp - nit000 + 1, cxios_context )   ! say to iom that we are at time step kstp 
     100                                CALL trc_rst_opn( istp )         ! Open tracer                                !   restart file 
     101                                CALL dta_dyn    ( istp )         ! Interpolation of the dynamical fields 
     102                                CALL trc_stp    ( istp )         ! time-stepping 
     103         IF( lk_vvl )           CALL dta_dyn_swp( istp )         ! swap of sea surface height and vertical scale factors 
     104                                CALL stp_ctl    ( istp, indic )  ! Time loop: control and print 
    105105         istp = istp + 1 
    106106         IF( lk_mpp )   CALL mpp_max( nstop ) 
     
    265265      IF( nn_timing == 1 )  CALL timing_start( 'nemo_init') 
    266266      ! 
    267                             CALL     phy_cst    ! Physical constants 
    268                             CALL     eos_init   ! Equation of state 
    269       IF( lk_c1d        )   CALL     c1d_init   ! 1D column configuration 
    270                             CALL     dom_cfg    ! Domain configuration 
     267                            CALL  phy_cst    ! Physical constants 
     268                            CALL  eos_init   ! Equation of state 
     269      IF( lk_c1d        )   CALL  c1d_init   ! 1D column configuration 
     270                            CALL  dom_cfg    ! Domain configuration 
     271      ! 
    271272      ! 
    272273      INQUIRE( FILE='coordinates.nc', EXIST = llexist )   ! Check if coordinate file exist 
    273274      ! 
    274       IF( llexist )  THEN  ;  CALL  dom_init   !  compute the grid from coordinates and bathymetry 
    275       ELSE                 ;  CALL  dom_rea    !  read grid from the meskmask 
     275      IF( llexist )  THEN ;  CALL  dom_init   !  compute the grid from coordinates and bathymetry 
     276      ELSE                ;  CALL  dom_rea    !  read grid from the meskmask 
    276277      ENDIF 
    277278                            CALL  istate_init   ! ocean initial state (Dynamics and tracers) 
    278279 
    279       IF( ln_nnogather )    CALL nemo_northcomms   ! Initialise the northfold neighbour lists (must be done after the masks are defined) 
    280  
    281       IF( ln_ctl        )   CALL prt_ctl_init   ! Print control 
     280      IF( ln_nnogather )    CALL  nemo_northcomms   ! Initialise the northfold neighbour lists (must be done after the masks are defined) 
     281 
     282      IF( ln_ctl       )    CALL prt_ctl_init   ! Print control 
    282283 
    283284                            CALL     sbc_init   ! Forcings : surface module 
     
    289290 
    290291                            CALL tra_qsr_init   ! penetrative solar radiation qsr 
    291       IF( lk_trabbl     )   CALL tra_bbl_init   ! advective (and/or diffusive) bottom boundary layer scheme 
     292      IF( lk_trabbl      CALL tra_bbl_init   ! advective (and/or diffusive) bottom boundary layer scheme 
    292293 
    293294                            CALL trc_nam_run    ! Needed to get restart parameters for passive tracers 
    294295                            CALL trc_rst_cal( nit000, 'READ' )   ! calendar 
    295296                            CALL dta_dyn_init   ! Initialization for the dynamics 
    296  
    297297                            CALL     trc_init   ! Passive tracers initialization 
    298                             CALL dia_ptr_init   ! Initialise diaptr as some variables are used  
    299298      !                                         ! in various advection and diffusion routines 
    300299      IF(lwp) WRITE(numout,cform_aaa)           ! Flag AAAAAAA 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r5363 r7522  
    2323   USE dom_oce         ! domain: ocean 
    2424   USE sbc_oce         ! surface boundary condition: ocean 
     25   USE trc_oce         ! shared ocean-passive tracers variables 
    2526   USE phycst          ! physical constants 
    2627   USE closea          ! closed seas 
     
    9798      END DO 
    9899      ! 
    99       IF( lk_vvl )           CALL dom_vvl_init ! Vertical variable mesh 
    100       ! 
    101       IF( lk_c1d         )   CALL cor_c1d      ! 1D configuration: Coriolis set at T-point 
    102       ! 
    103       ! 
    104       hu(:,:) = 0._wp                          ! Ocean depth at U-points 
    105       hv(:,:) = 0._wp                          ! Ocean depth at V-points 
    106       ht(:,:) = 0._wp                          ! Ocean depth at T-points 
    107       DO jk = 1, jpkm1 
    108          hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk) 
    109          hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk) 
    110          ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 
    111       END DO 
    112       !                                        ! Inverse of the local depth 
    113       hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:) 
    114       hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:) 
     100      IF( lk_c1d )           CALL cor_c1d      ! 1D configuration: Coriolis set at T-point 
     101      ! 
     102      IF( .NOT.lk_offline ) THEN 
     103        ! 
     104        IF( lk_vvl )         CALL dom_vvl_init ! Vertical variable mesh 
     105        ! 
     106        hu(:,:) = 0._wp                          ! Ocean depth at U-points 
     107        hv(:,:) = 0._wp                          ! Ocean depth at V-points 
     108        ht(:,:) = 0._wp                          ! Ocean depth at T-points 
     109        DO jk = 1, jpkm1 
     110           hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk) 
     111           hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk) 
     112           ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 
     113        END DO 
     114        !                                        ! Inverse of the local depth 
     115        hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:) 
     116        hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:) 
     117        ! 
     118      ENDIF 
    115119 
    116120                             CALL dom_stp      ! time step 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90

    r6969 r7522  
    263263      ENDIF 
    264264 
    265       ! set the number of level over which river runoffs are applied  
    266       ! online configuration : computed in sbcrnf 
    267       IF( lk_offline ) THEN 
    268         nk_rnf(:,:) = 1 
    269         h_rnf (:,:) = e3t_0(:,:,1) 
    270       ENDIF 
    271265 
    272266      ! dust input from the atmosphere 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/TRP/trcadv.F90

    r5385 r7522  
    8888         r2dt(:) = 2. * rdttrc(:)       ! = 2 rdttrc (leapfrog) 
    8989      ENDIF 
    90       !                                                   ! effective transport 
    91       DO jk = 1, jpkm1 
    92          !                                                ! eulerian transport only 
    93          zun(:,:,jk) = e2u  (:,:) * fse3u(:,:,jk) * un(:,:,jk) 
    94          zvn(:,:,jk) = e1v  (:,:) * fse3v(:,:,jk) * vn(:,:,jk) 
    95          zwn(:,:,jk) = e1e2t(:,:)                 * wn(:,:,jk) 
    96          ! 
    97       END DO 
    98       ! 
    99       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 
    100          zun(:,:,:) = zun(:,:,:) + un_td(:,:,:) 
    101          zvn(:,:,:) = zvn(:,:,:) + vn_td(:,:,:) 
     90      !   
     91      IF( lk_offline ) THEN 
     92         zun(:,:,:) = un(:,:,:)     ! effective transport already in un/vn/wn 
     93         zvn(:,:,:) = vn(:,:,:) 
     94         zwn(:,:,:) = wn(:,:,:) 
     95      ELSE 
     96         !                                                         ! effective transport 
     97         DO jk = 1, jpkm1 
     98            !                                                ! eulerian transport only 
     99            zun(:,:,jk) = e2u  (:,:) * fse3u(:,:,jk) * un(:,:,jk) 
     100            zvn(:,:,jk) = e1v  (:,:) * fse3v(:,:,jk) * vn(:,:,jk) 
     101            zwn(:,:,jk) = e1e2t(:,:)                 * wn(:,:,jk) 
     102            ! 
     103         END DO 
     104         ! 
     105         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 
     106            zun(:,:,:) = zun(:,:,:) + un_td(:,:,:) 
     107            zvn(:,:,:) = zvn(:,:,:) + vn_td(:,:,:) 
     108         ENDIF 
     109         ! 
     110         zun(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom 
     111         zvn(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom 
     112         zwn(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom 
     113         ! 
     114 
     115         IF( lk_traldf_eiv .AND. .NOT. ln_traldf_grif )   &  ! add the eiv transport (if necessary) 
     116            &              CALL tra_adv_eiv( kt, nittrc000, zun, zvn, zwn, 'TRC' ) 
     117         ! 
     118         IF( ln_mle    )   CALL tra_adv_mle( kt, nittrc000, zun, zvn, zwn, 'TRC' )    ! add the mle transport (if necessary) 
     119         ! 
    102120      ENDIF 
    103       ! 
    104       zun(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom 
    105       zvn(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom 
    106       zwn(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom 
    107  
    108       IF( lk_traldf_eiv .AND. .NOT. ln_traldf_grif )   &  ! add the eiv transport (if necessary) 
    109          &              CALL tra_adv_eiv( kt, nittrc000, zun, zvn, zwn, 'TRC' ) 
    110       ! 
    111       IF( ln_mle    )   CALL tra_adv_mle( kt, nittrc000, zun, zvn, zwn, 'TRC' )    ! add the mle transport (if necessary) 
    112121      ! 
    113122      SELECT CASE ( nadv )                            !==  compute advection trend and add it to general trend  ==! 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/TRP/trcrad.F90

    r7494 r7522  
    6565      IF( lk_c14b    )   CALL trc_rad_sms( kt, trb, trn, jp_c14b0, jp_c14b1              )  ! bomb C14 
    6666      IF( lk_pisces  )   CALL trc_rad_sms( kt, trb, trn, jp_pcs0 , jp_pcs1, cpreserv='Y' )  ! PISCES model 
    67       IF( lk_my_trc  )   CALL trc_rad_sms( kt, trb, trn, jp_myt0 , jp_myt1              )  ! MY_TRC model 
     67      IF( lk_my_trc  )   CALL trc_rad_sms( kt, trb, trn, jp_myt0 , jp_myt1, cpreserv='Y' )  ! MY_TRC model 
    6868 
    6969      ! 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90

    r6941 r7522  
    128128      ! Coupling offline : runoff are in emp which contains E-P-R 
    129129      ! 
    130       IF( .NOT. lk_offline .AND. lk_vvl ) THEN  ! online coupling with vvl 
     130      IF( lk_vvl ) THEN                         ! linear free surface vvl 
    131131         zsfx(:,:) = 0._wp 
    132       ELSE                                      ! online coupling free surface or offline with free surface 
     132      ELSE                                      ! no vvl 
    133133         zsfx(:,:) = emp(:,:) 
    134134      ENDIF 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/trcstp.F90

    r6967 r7522  
    8686         tra(:,:,:,:) = 0.e0 
    8787         ! 
    88                                    CALL trc_rst_opn  ( kt )       ! Open tracer restart file  
     88         IF( .NOT.lk_offline )     CALL trc_rst_opn  ( kt )       ! Open tracer restart file  
    8989         IF( lrst_trc )            CALL trc_rst_cal  ( kt, 'WRITE' )   ! calendar 
    9090         IF( lk_iomput ) THEN  ;   CALL trc_wri      ( kt )       ! output of passive tracers with iom I/O manager 
Note: See TracChangeset for help on using the changeset viewer.