New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 6952 – NEMO

Changeset 6952


Ignore:
Timestamp:
2016-09-23T16:02:11+02:00 (8 years ago)
Author:
cetlod
Message:

Offline vvl: 1st implementation, see ticket #1775

Location:
branches/2015/dev_r6946_Offline_vvl/NEMOGCM
Files:
13 added
11 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/OFF_SRC/domrea.F90

    r5504 r6952  
    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/dev_r6946_Offline_vvl/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90

    r5767 r6952  
    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 
    63    INTEGER  , SAVE      ::   jf_emp         ! index of water flux 
     68   INTEGER  , SAVE      ::   jf_empn        ! 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 
    9287 
    9388   !! * Substitutions 
     
    113108      !!---------------------------------------------------------------------- 
    114109      ! 
    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       ! 
     110      USE oce, ONLY:  zhdivtr => ua 
    120111      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    !  
     112      INTEGER             ::   ji, jj, jk 
     113      REAL(wp), POINTER, DIMENSION(:,:)   :: zemp 
     114      ! 
    127115      !!---------------------------------------------------------------------- 
    128116       
     
    130118      IF( nn_timing == 1 )  CALL timing_start( 'dta_dyn') 
    131119      ! 
    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 
     120      CALL fld_read( kt, 1, sf_dyn )      !=  read data at kt time step   ==! 
    243121      ! 
    244122      tsn(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:)  * tmask(:,:,:)    ! temperature 
    245123      tsn(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:)  * tmask(:,:,:)    ! salinity 
    246       ! 
     124      wndm(:,:)         = sf_dyn(jf_wnd)%fnow(:,:,1)  * tmask(:,:,1)    ! wind speed - needed for gas exchange 
     125      fmmflx(:,:)       = sf_dyn(jf_fmf)%fnow(:,:,1)  * tmask(:,:,1)    ! downward salt flux (v3.5+) 
     126      fr_i(:,:)         = sf_dyn(jf_ice)%fnow(:,:,1)  * tmask(:,:,1)    ! Sea-ice fraction 
     127      qsr (:,:)         = sf_dyn(jf_qsr)%fnow(:,:,1)  * tmask(:,:,1)    ! solar radiation 
     128      emp (:,:)         = sf_dyn(jf_empn)%fnow(:,:,1) * tmask(:,:,1)    ! E-P 
     129      emp_b (:,:)       = sf_dyn(jf_empb)%fnow(:,:,1) * tmask(:,:,1)    ! E-P 
     130      IF( ln_dynrnf ) THEN  
     131         rnf (:,:)      = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask(:,:,1)    ! E-P 
     132         IF( ln_dynrnf_depth .AND. lk_vvl )    CALL  dta_dyn_hrnf 
     133      ENDIF 
     134      ! 
     135      un(:,:,:)        = sf_dyn(jf_uwd)%fnow(:,:,:) * umask(:,:,:)    ! effective u-transport 
     136      vn(:,:,:)        = sf_dyn(jf_vwd)%fnow(:,:,:) * vmask(:,:,:)    ! effective v-transport 
     137      wn(:,:,:)        = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:)    ! effective v-transport 
     138      zhdivtr(:,:,:)   = sf_dyn(jf_div)%fnow(:,:,:) * tmask(:,:,:)    ! effective u-transport 
     139      ! 
     140      IF( lk_vvl ) THEN 
     141         CALL wrk_alloc(jpi, jpj, zemp ) 
     142         zemp(:,:) = 0.5_wp * ( emp(:,:) + emp_b(:,:) ) + rnf(:,:) + fwbcorr * tmask(:,:,1) 
     143         CALL dta_dyn_ssh( kt, zhdivtr, sshb, zemp, ssha, fse3t_a(:,:,:) )  !=  ssh, vertical scale factor & vertical transport 
     144         CALL wrk_dealloc(jpi, jpj, zemp ) 
     145         !                                           Write in the tracer restart file 
     146         !                                          ******************************* 
     147         IF( lrst_trc ) THEN 
     148            IF(lwp) WRITE(numout,*) 
     149            IF(lwp) WRITE(numout,*) 'dta_dyn_ssh : ssh field written in tracer restart file ',   & 
     150               &                    'at it= ', kt,' date= ', ndastp 
     151            IF(lwp) WRITE(numout,*) '~~~~' 
     152            CALL iom_rstput( kt, nitrst, numrtw, 'sshn', ssha ) 
     153            CALL iom_rstput( kt, nitrst, numrtw, 'sshb', sshn ) 
     154         ENDIF 
     155      ENDIF 
    247156      ! 
    248157      CALL eos    ( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 
     
    251160 
    252161      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 
     162      CALL zdf_mxl( kt )                                                   ! In any case, we need mxl 
     163      ! 
     164      IF( lk_ldfslp .AND. .NOT.lk_c1d )   CALL  dta_dyn_slp( kt )    ! Computation of slopes 
     165      ! 
     166      hmld(:,:)         = sf_dyn(jf_mld)%fnow(:,:,1) * tmask(:,:,1)    ! mixed layer depht 
     167      avt(:,:,:)        = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:)    ! vertical diffusive coefficient  
     168      ! 
    270169#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 
     170      ahu_bbl(:,:)      = sf_dyn(jf_ubl)%fnow(:,:,1) * umask(:,:,1)    ! bbl diffusive coef 
     171      ahv_bbl(:,:)      = sf_dyn(jf_vbl)%fnow(:,:,1) * vmask(:,:,1) 
    278172#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 
     173      ! 
     174      ! 
     175      CALL eos( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 
    302176      ! 
    303177      IF(ln_ctl) THEN                  ! print control 
     
    308182         CALL prt_ctl(tab3d_1=wn               , clinfo1=' wn      - : ', mask1=tmask, ovlap=1, kdim=jpk   ) 
    309183         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 ) 
     184!         CALL prt_ctl(tab2d_1=fr_i             , clinfo1=' fr_i    - : ', mask1=tmask, ovlap=1 ) 
     185!         CALL prt_ctl(tab2d_1=hmld             , clinfo1=' hmld    - : ', mask1=tmask, ovlap=1 ) 
     186!         CALL prt_ctl(tab2d_1=fmmflx           , clinfo1=' fmmflx  - : ', mask1=tmask, ovlap=1 ) 
     187!         CALL prt_ctl(tab2d_1=emp              , clinfo1=' emp     - : ', mask1=tmask, ovlap=1 ) 
     188!         CALL prt_ctl(tab2d_1=wndm             , clinfo1=' wspd    - : ', mask1=tmask, ovlap=1 ) 
     189!         CALL prt_ctl(tab2d_1=qsr              , clinfo1=' qsr     - : ', mask1=tmask, ovlap=1 ) 
    316190      ENDIF 
    317191      ! 
     
    335209      INTEGER  :: inum, idv, idimv                   ! local integer 
    336210      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 
     211      INTEGER  :: ji, jj, jk 
     212      REAL(wp) :: zcoef 
     213      INTEGER  :: nkrnf_max 
     214      REAL(wp) :: hrnf_max 
     215      !! 
     216      CHARACTER(len=100)            ::  cn_dir        !   Root directory for location of core files 
     217      TYPE(FLD_N), DIMENSION(jpfld) ::  slf_d         ! array of namelist informations on the fields to read 
     218      TYPE(FLD_N) :: sn_uwd, sn_vwd, sn_wwd, sn_empb, sn_empn  ! informations about the fields to be read 
     219      TYPE(FLD_N) :: sn_tem , sn_sal , sn_avt   !   "                 " 
     220      TYPE(FLD_N) :: sn_mld, sn_qsr, sn_wnd , sn_ice , sn_fmf   !   "               " 
     221      TYPE(FLD_N) :: sn_ubl, sn_vbl, sn_rnf    !   "              " 
     222      TYPE(FLD_N) :: sn_div  ! informations about the fields to be read 
     223 
     224      !!---------------------------------------------------------------------- 
     225      ! 
     226      NAMELIST/namdta_dyn/cn_dir, ln_dynrnf, ln_dynrnf_depth,  fwbcorr, & 
     227         &                sn_uwd, sn_vwd, sn_wwd, sn_empb, sn_empn,   & 
     228         &                sn_avt, sn_tem, sn_sal, sn_mld , sn_qsr ,   & 
     229         &                sn_wnd, sn_ice, sn_fmf,                    & 
     230         &                sn_ubl, sn_vbl, sn_rnf,                   & 
     231         &                sn_div  
    349232      ! 
    350233      REWIND( numnam_ref )              ! Namelist namdta_dyn in reference namelist : Offline: init. of dynamical data 
     
    363246         WRITE(numout,*) '~~~~~~~ ' 
    364247         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 
     248         WRITE(numout,*) '      runoffs option enabled (T) or not (F)            ln_dynrnf        = ', ln_dynrnf 
     249         WRITE(numout,*) '      runoffs is spread in vertical                    ln_dynrnf_depth  = ', ln_dynrnf_depth 
     250         WRITE(numout,*) '      annual global mean of empmr for ssh correction   fwbcorr          = ', fwbcorr 
    369251         WRITE(numout,*) 
    370252      ENDIF 
    371253      !  
    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  
     254 
     255      jf_uwd  = 1     ;   jf_vwd  = 2    ;   jf_wwd = 3   
     256      jf_empb = 4     ;   jf_empn = 5    ;   jf_avt = 6 
     257      jf_tem  = 7     ;   jf_sal  = 8    ;   jf_mld = 9    ;   jf_qsr = 10 
     258      jf_wnd  = 11    ;   jf_ice  = 12   ;   jf_fmf = 13   ;   jfld   = jf_fmf 
     259 
     260      ! 
     261      slf_d(jf_uwd)  = sn_uwd    ;   slf_d(jf_vwd)  = sn_vwd   ;   slf_d(jf_wwd) = sn_wwd 
     262      slf_d(jf_empb) = sn_empb   ;   slf_d(jf_empn) = sn_empn  ;   slf_d(jf_avt) = sn_avt 
     263      slf_d(jf_tem)  = sn_tem    ;   slf_d(jf_sal)  = sn_sal   ;   slf_d(jf_mld) = sn_mld 
     264      slf_d(jf_qsr)  = sn_qsr    ;   slf_d(jf_wnd)  = sn_wnd   ;   slf_d(jf_ice) = sn_ice 
     265      slf_d(jf_fmf)  = sn_fmf 
     266 
     267      ! 
     268      IF( lk_vvl ) THEN 
     269                 jf_div  = jfld + 1    ;     jfld = jf_div 
     270           slf_d(jf_div) = sn_div     
     271      ENDIF 
     272      ! 
     273      IF( lk_trabbl ) THEN 
     274                 jf_ubl  = jfld + 1    ;         jf_vbl  = jfld + 2     ;      jfld = jf_vbl 
     275           slf_d(jf_ubl) = sn_ubl      ;   slf_d(jf_vbl) = sn_vbl 
     276      ENDIF 
    389277      ! 
    390278      IF( ln_dynrnf ) THEN 
    391                 jf_rnf = jfld + 1  ;  jfld  = jf_rnf 
    392          slf_d(jf_rnf) = sn_rnf 
     279                jf_rnf  = jfld + 1     ;     jfld  = jf_rnf 
     280          slf_d(jf_rnf) = sn_rnf     
    393281      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 
     282         rnf(:,:) = 0._wp 
     283      ENDIF 
     284 
    428285   
    429286      ALLOCATE( sf_dyn(jfld), STAT=ierr )         ! set sf structure 
    430       IF( ierr > 0 ) THEN 
     287      IF( ierr > 0 )  THEN 
    431288         CALL ctl_stop( 'dta_dyn: unable to allocate sf structure' )   ;   RETURN 
    432289      ENDIF 
    433290      !                                         ! fill sf with slf_i and control print 
    434291      CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' ) 
     292      ! 
    435293      ! Open file for each variable to get his number of dimension 
    436294      DO ifpr = 1, jfld 
     
    456314            ALLOCATE( uslpdta (jpi,jpj,jpk,2), vslpdta (jpi,jpj,jpk,2),    & 
    457315            &         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 
     316            ! 
     317            IF( ierr2 > 0 )  THEN 
     318               CALL ctl_stop( 'dta_dyn_init : unable to allocate slope arrays' )   ;   RETURN 
     319            ENDIF 
    464320         ENDIF 
    465321      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  
     322      ! 
     323      IF( lk_vvl ) THEN 
     324        IF( .NOT. sf_dyn(jf_uwd)%ln_clim .AND. ln_rsttr .AND.    &                     ! Restart: read in restart file 
     325           iom_varid( numrtr, 'sshn', ldstop = .FALSE. ) > 0 ) THEN 
     326           IF(lwp) WRITE(numout,*) ' sshn forcing fields read in the restart file for initialisation' 
     327           CALL iom_get( numrtr, jpdom_autoglo, 'sshn', sshn(:,:)   ) 
     328           CALL iom_get( numrtr, jpdom_autoglo, 'sshb', sshb(:,:)   ) 
     329        ELSE 
     330           IF(lwp) WRITE(numout,*) ' sshn forcing fields read in the restart file for initialisation' 
     331           CALL iom_open( 'restart', inum ) 
     332           CALL iom_get( inum, jpdom_autoglo, 'sshn', sshn(:,:)   ) 
     333           CALL iom_get( inum, jpdom_autoglo, 'sshb', sshb(:,:)   ) 
     334           CALL iom_close( inum )                                        ! close file 
     335        ENDIF 
     336        ! 
     337        DO jk = 1, jpkm1 
     338           fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) ) 
     339        ENDDO 
     340        fse3t_a(:,:,jpk) = e3t_0(:,:,jpk) 
     341 
     342        ! Horizontal scale factor interpolations 
     343        ! -------------------------------------- 
     344        CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' ) 
     345        CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 
     346 
     347        ! Vertical scale factor interpolations 
     348        ! ------------------------------------ 
     349        CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n(:,:,:), 'W' ) 
     350   
     351        fse3t_b(:,:,:)  = fse3t_n(:,:,:) 
     352        fse3u_b(:,:,:)  = fse3u_n(:,:,:) 
     353        fse3v_b(:,:,:)  = fse3v_n(:,:,:) 
     354 
     355        ! t- and w- points depth 
     356        ! ---------------------- 
     357        fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 
     358        fsdepw_n(:,:,1) = 0.0_wp 
     359 
     360        DO jk = 2, jpk 
     361           DO jj = 1,jpj 
     362              DO ji = 1,jpi 
     363                !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere 
     364                !    tmask = wmask, ie everywhere expect at jk = mikt 
     365                                                                   ! 1 for jk = 
     366                                                                   ! mikt 
     367                 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
     368                 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 
     369                 fsdept_n(ji,jj,jk) =      zcoef  * ( fsdepw_n(ji,jj,jk  ) + 0.5 * fse3w_n(ji,jj,jk))  & 
     370                     &                + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)) 
     371              END DO 
     372           END DO 
     373        END DO 
     374 
     375        fsdept_b(:,:,:) = fsdept_n(:,:,:) 
     376        fsdepw_b(:,:,:) = fsdepw_n(:,:,:) 
     377        ! 
     378      ENDIF 
     379      ! 
     380      IF( ln_dynrnf .AND. ln_dynrnf_depth ) THEN       ! read depht over which runoffs are distributed 
     381         IF(lwp) WRITE(numout,*)  
     382         IF(lwp) WRITE(numout,*) ' read in the file depht over which runoffs are distributed' 
     383         CALL iom_open ( "runoffs", inum )                           ! open file 
     384         CALL iom_get  ( inum, jpdom_data, 'rodepth', h_rnf )   ! read the river mouth array 
     385         CALL iom_close( inum )                                        ! close file 
     386         ! 
     387         nk_rnf(:,:) = 0                               ! set the number of level over which river runoffs are applied 
     388         DO jj = 1, jpj 
     389            DO ji = 1, jpi 
     390               IF( h_rnf(ji,jj) > 0._wp ) THEN 
     391                  jk = 2 
     392                  DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1 
     393                  END DO 
     394                  nk_rnf(ji,jj) = jk 
     395               ELSEIF( h_rnf(ji,jj) == -1._wp   ) THEN   ;  nk_rnf(ji,jj) = 1 
     396               ELSEIF( h_rnf(ji,jj) == -999._wp ) THEN   ;  nk_rnf(ji,jj) = mbkt(ji,jj) 
     397               ELSE 
     398                  CALL ctl_stop( 'sbc_rnf_init: runoff depth not positive, and not -999 or -1, rnf value in file fort.999'  ) 
     399                  WRITE(999,*) 'ji, jj, h_rnf(ji,jj) :', ji, jj, h_rnf(ji,jj) 
     400               ENDIF 
    515401            END DO 
    516402         END DO 
     403         DO jj = 1, jpj                                ! set the associated depth 
     404            DO ji = 1, jpi 
     405               h_rnf(ji,jj) = 0._wp 
     406               DO jk = 1, nk_rnf(ji,jj) 
     407                  h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk) 
     408               END DO 
     409            END DO 
     410         END DO 
     411      ELSE                                       ! runoffs applied at the surface 
     412         nk_rnf(:,:) = 1 
     413         h_rnf (:,:) = fse3t(:,:,1) 
     414      ENDIF 
     415      nkrnf_max = MAXVAL( nk_rnf(:,:) ) 
     416      hrnf_max = MAXVAL( h_rnf(:,:) ) 
     417      IF( lk_mpp )  THEN 
     418         CALL mpp_max( nkrnf_max )                 ! max over the  global domain 
     419         CALL mpp_max( hrnf_max )                 ! max over the  global domain 
     420      ENDIF 
     421      IF(lwp) WRITE(numout,*) ' ' 
     422      IF(lwp) WRITE(numout,*) ' max depht of runoff : ', hrnf_max,'    max level  : ', nkrnf_max 
     423      IF(lwp) WRITE(numout,*) ' ' 
     424      ! 
     425      CALL dta_dyn( nit000 ) 
     426      ! 
     427   END SUBROUTINE dta_dyn_init 
     428 
     429   SUBROUTINE dta_dyn_swp( kt ) 
     430     !!--------------------------------------------------------------------- 
     431      !!                    ***  ROUTINE dta_dyn_swp  *** 
     432      !! 
     433      !! ** Purpose : Swap and the data and compute the vertical scale factor at U/V/W point 
     434      !!              and the depht 
     435      !! 
     436      !!--------------------------------------------------------------------- 
     437      INTEGER, INTENT(in) :: kt       ! time step 
     438      INTEGER             :: ji, jj, jk 
     439      REAL(wp)            :: zcoef 
     440      ! 
     441      !!--------------------------------------------------------------------- 
     442 
     443      IF( kt == nit000 ) THEN 
     444         IF(lwp) WRITE(numout,*) 
     445         IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height' 
     446         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
     447      ENDIF 
     448 
     449      sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:))  ! before <-- now filtered 
     450      sshn(:,:) = ssha(:,:) 
     451 
     452      fse3t_n(:,:,:) = fse3t_a(:,:,:) 
     453 
     454      ! Reconstruction of all vertical scale factors at now and before time steps 
     455      ! ============================================================================= 
     456 
     457      ! Horizontal scale factor interpolations 
     458      ! -------------------------------------- 
     459      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' ) 
     460      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 
     461 
     462      ! Vertical scale factor interpolations 
     463      ! ------------------------------------ 
     464      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' ) 
     465 
     466      fse3t_b(:,:,:)  = fse3t_n(:,:,:) 
     467      fse3u_b(:,:,:)  = fse3u_n(:,:,:) 
     468      fse3v_b(:,:,:)  = fse3v_n(:,:,:) 
     469 
     470      ! t- and w- points depth 
     471      ! ---------------------- 
     472      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 
     473      fsdepw_n(:,:,1) = 0.0_wp 
     474 
     475      DO jk = 2, jpk 
     476         DO jj = 1,jpj 
     477            DO ji = 1,jpi 
     478                 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
     479                 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 
     480                 fsdept_n(ji,jj,jk) =      zcoef  * ( fsdepw_n(ji,jj,jk  ) + 0.5 * fse3w_n(ji,jj,jk))  & 
     481                     &                + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)) 
     482              END DO 
     483           END DO 
     484        END DO 
     485 
     486      fsdept_b(:,:,:) = fsdept_n(:,:,:) 
     487      fsdepw_b(:,:,:) = fsdepw_n(:,:,:) 
     488 
     489      ! 
     490   END SUBROUTINE dta_dyn_swp 
     491 
     492   SUBROUTINE dta_dyn_ssh( kt, phdivtr, psshb,  pemp, pssha, pe3ta ) 
     493      !!---------------------------------------------------------------------- 
     494      !!                ***  ROUTINE dta_dyn_wzv  *** 
     495      !!                    
     496      !! ** Purpose :   compute the after ssh (ssha) and the now vertical velocity 
     497      !! 
     498      !! ** Method  : Using the incompressibility hypothesis,  
     499      !!        - the ssh increment is computed by integrating the horizontal divergence  
     500      !!          and multiply by the time step. 
     501      !! 
     502      !!        - compute the after scale factor : repartition of ssh INCREMENT proportionnaly 
     503      !!                                           to the level thickness ( z-star case ) 
     504      !! 
     505      !!        - the vertical velocity is computed by integrating the horizontal divergence   
     506      !!          from the bottom to the surface minus the scale factor evolution. 
     507      !!          The boundary conditions are w=0 at the bottom (no flux) 
     508      !! 
     509      !! ** action  :   ssha / e3t_a / wn 
     510      !! 
     511      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
     512      !!---------------------------------------------------------------------- 
     513      !! * Arguments 
     514      INTEGER,                                   INTENT(in )    :: kt        !  time-step 
     515      REAL(wp), DIMENSION(jpi,jpj,jpk)          , INTENT(in )   :: phdivtr   ! horizontal divergence transport 
     516      REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(in )   :: psshb     ! now ssh 
     517      REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(in )   :: pemp      ! evaporation minus precipitation 
     518      REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(inout) :: pssha     ! after ssh 
     519      REAL(wp), DIMENSION(jpi,jpj,jpk), OPTIONAL, INTENT(out)   :: pe3ta     ! after vertical scale factor 
     520      !! * Local declarations 
     521      INTEGER                       :: jk 
     522      REAL(wp), DIMENSION(jpi,jpj)  :: zhdiv   
     523      REAL(wp)  :: z2dt   
     524      !!---------------------------------------------------------------------- 
     525       
     526      ! 
     527      z2dt = 2._wp * rdt 
     528      ! 
     529      zhdiv(:,:) = 0._wp 
     530      DO jk = 1, jpkm1 
     531         zhdiv(:,:) = zhdiv(:,:) +  phdivtr(:,:,jk) * tmask(:,:,jk) 
    517532      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) 
     533      !                                                ! Sea surface  elevation time-stepping 
     534      pssha(:,:) = ( psshb(:,:) - z2dt * ( r1_rau0 * pemp(:,:)  + zhdiv(:,:) ) ) * ssmask(:,:) 
     535      !                                                 !  
     536      !                                                 ! After acale factors at t-points ( z_star coordinate ) 
     537      DO jk = 1, jpkm1 
     538        pe3ta(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + pssha(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) ) 
    524539      END DO 
    525540      ! 
    526    END SUBROUTINE dta_dyn_wzv 
    527  
    528    SUBROUTINE dta_dyn_slp( kt, pts, puslp, pvslp, pwslpi, pwslpj ) 
     541   END SUBROUTINE dta_dyn_ssh 
     542 
     543 
     544   SUBROUTINE dta_dyn_hrnf 
     545      !!---------------------------------------------------------------------- 
     546      !!                  ***  ROUTINE sbc_rnf  *** 
     547      !! 
     548      !! ** Purpose :   update the horizontal divergence with the runoff inflow 
     549      !! 
     550      !! ** Method  : 
     551      !!                CAUTION : rnf is positive (inflow) decreasing the 
     552      !!                          divergence and expressed in m/s 
     553      !! 
     554      !! ** Action  :   phdivn   decreased by the runoff inflow 
     555      !!---------------------------------------------------------------------- 
     556      !! 
     557      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     558      !!---------------------------------------------------------------------- 
     559      ! 
     560      DO jj = 1, jpj                   ! update the depth over which runoffs are distributed 
     561         DO ji = 1, jpi 
     562            h_rnf(ji,jj) = 0._wp 
     563            DO jk = 1, nk_rnf(ji,jj)                           ! recalculates h_rnf to be the depth in metres 
     564                h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk)   ! to the bottom of the relevant grid box 
     565            END DO 
     566        END DO 
     567      END DO 
     568      ! 
     569   END SUBROUTINE dta_dyn_hrnf 
     570 
     571 
     572 
     573   SUBROUTINE dta_dyn_slp( kt ) 
     574      !!--------------------------------------------------------------------- 
     575      !!                    ***  ROUTINE dta_dyn_slp  *** 
     576      !! 
     577      !! ** Purpose : Computation of slope 
     578      !! 
     579      !!--------------------------------------------------------------------- 
     580      USE oce, ONLY:  zts => tsa  
     581      ! 
     582      INTEGER,  INTENT(in) :: kt       ! time step 
     583      ! 
     584      INTEGER  ::   ji, jj     ! dummy loop indices 
     585      REAL(wp) ::   ztinta     ! ratio applied to after  records when doing time interpolation 
     586      REAL(wp) ::   ztintb     ! ratio applied to before records when doing time interpolation 
     587      INTEGER  ::   iswap, isecdyn, iprevrec  
     588      REAL(wp), POINTER, DIMENSION(:,:,:) :: zuslp, zvslp, zwslpi, zwslpj 
     589      !!--------------------------------------------------------------------- 
     590      ! 
     591      CALL wrk_alloc(jpi, jpj, jpk, zuslp, zvslp, zwslpi, zwslpj ) 
     592      ! 
     593      isecdyn = nsec_year + nsec1jan000   ! number of seconds between Jan. 1st 00h of nit000 year and the middle of time step 
     594      ! 
     595      IF( kt == nit000 ) THEN    ;    iprevrec = 0 
     596      ELSE                       ;    iprevrec = sf_dyn(jf_tem)%nrec_a(2) 
     597      ENDIF 
     598      ! 
     599      IF( sf_dyn(jf_tem)%ln_tint ) THEN    ! Computes slopes (here avt is used as workspace)                        
     600         IF( kt == nit000 ) THEN 
     601            zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,1) * tmask(:,:,:)   ! temperature 
     602            zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:)   ! salinity  
     603            avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,1) * tmask(:,:,:)   ! vertical diffusive coef. 
     604            CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
     605            uslpdta (:,:,:,1) = zuslp (:,:,:)  
     606            vslpdta (:,:,:,1) = zvslp (:,:,:)  
     607            wslpidta(:,:,:,1) = zwslpi(:,:,:)  
     608            wslpjdta(:,:,:,1) = zwslpj(:,:,:)  
     609            ! 
     610            zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:)   ! temperature 
     611            zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity  
     612            avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef. 
     613            CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
     614            uslpdta (:,:,:,2) = zuslp (:,:,:)  
     615            vslpdta (:,:,:,2) = zvslp (:,:,:)  
     616            wslpidta(:,:,:,2) = zwslpi(:,:,:)  
     617            wslpjdta(:,:,:,2) = zwslpj(:,:,:)  
     618         ELSE 
     619           !  
     620           iswap = 0 
     621           IF( sf_dyn(jf_tem)%nrec_a(2) - iprevrec /= 0 )  iswap = 1 
     622           IF( isecdyn > sf_dyn(jf_tem)%nrec_b(2) .AND. iswap == 1 )  THEN    ! read/update the after data 
     623              IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt 
     624              uslpdta (:,:,:,1) =  uslpdta (:,:,:,2)         ! swap the data 
     625              vslpdta (:,:,:,1) =  vslpdta (:,:,:,2)   
     626              wslpidta(:,:,:,1) =  wslpidta(:,:,:,2)  
     627              wslpjdta(:,:,:,1) =  wslpjdta(:,:,:,2)  
     628              ! 
     629              zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:)   ! temperature 
     630              zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity  
     631              avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef. 
     632              CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
     633              ! 
     634              uslpdta (:,:,:,2) = zuslp (:,:,:)  
     635              vslpdta (:,:,:,2) = zvslp (:,:,:)  
     636              wslpidta(:,:,:,2) = zwslpi(:,:,:)  
     637              wslpjdta(:,:,:,2) = zwslpj(:,:,:)  
     638            ENDIF 
     639         ENDIF 
     640      ENDIF 
     641      ! 
     642      IF( sf_dyn(jf_tem)%ln_tint )  THEN 
     643         ztinta =  REAL( isecdyn - sf_dyn(jf_tem)%nrec_b(2), wp )  & 
     644            &    / REAL( sf_dyn(jf_tem)%nrec_a(2) - sf_dyn(jf_tem)%nrec_b(2), wp ) 
     645         ztintb =  1. - ztinta 
     646#if defined key_ldfslp && ! defined key_c1d 
     647         uslp (:,:,:) = ztintb * uslpdta (:,:,:,1)  + ztinta * uslpdta (:,:,:,2)   
     648         vslp (:,:,:) = ztintb * vslpdta (:,:,:,1)  + ztinta * vslpdta (:,:,:,2)   
     649         wslpi(:,:,:) = ztintb * wslpidta(:,:,:,1)  + ztinta * wslpidta(:,:,:,2)   
     650         wslpj(:,:,:) = ztintb * wslpjdta(:,:,:,1)  + ztinta * wslpjdta(:,:,:,2)   
     651#endif 
     652      ELSE 
     653         zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:)   ! temperature 
     654         zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:)   ! salinity  
     655         avt(:,:,:)        = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:)   ! vertical diffusive coef. 
     656         CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
     657         ! 
     658#if defined key_ldfslp && ! defined key_c1d 
     659         uslp (:,:,:) = zuslp (:,:,:) 
     660         vslp (:,:,:) = zvslp (:,:,:) 
     661         wslpi(:,:,:) = zwslpi(:,:,:) 
     662         wslpj(:,:,:) = zwslpj(:,:,:) 
     663#endif 
     664      ENDIF 
     665      ! 
     666      CALL wrk_dealloc(jpi, jpj, jpk, zuslp, zvslp, zwslpi, zwslpj ) 
     667      ! 
     668   END SUBROUTINE dta_dyn_slp 
     669 
     670   SUBROUTINE compute_slopes( kt, pts, puslp, pvslp, pwslpi, pwslpj ) 
    529671      !!--------------------------------------------------------------------- 
    530672      !!                    ***  ROUTINE dta_dyn_slp  *** 
     
    568710#endif 
    569711      ! 
    570    END SUBROUTINE dta_dyn_slp 
     712   END SUBROUTINE compute_slopes 
     713 
    571714   !!====================================================================== 
    572715END MODULE dtadyn 
  • branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/OFF_SRC/nemogcm.F90

    r5504 r6952  
    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/dev_r6946_Offline_vvl/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r5363 r6952  
    2323   USE dom_oce         ! domain: ocean 
    2424   USE sbc_oce         ! surface boundary condition: ocean 
     25   USE trc_oce, ONLY : lk_offline         !  
    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(:,:) 
    115  
     100      IF( .NOT. lk_offline ) THEN 
     101      ! 
     102         IF( lk_vvl         )      CALL dom_vvl_init ! Vertical variable mesh 
     103         ! 
     104         IF( lk_c1d         )      CALL cor_c1d      ! 1D configuration: Coriolis set at T-point 
     105         ! 
     106         ! 
     107         hu(:,:) = 0._wp                          ! Ocean depth at U-points 
     108         hv(:,:) = 0._wp                          ! Ocean depth at V-points 
     109         ht(:,:) = 0._wp                          ! Ocean depth at T-points 
     110         DO jk = 1, jpkm1 
     111            hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk) 
     112            hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk) 
     113            ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 
     114         END DO 
     115         !                                        ! Inverse of the local depth 
     116         hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:) 
     117         hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:) 
     118        ! 
     119      ENDIF 
    116120                             CALL dom_stp      ! time step 
    117121      IF( nmsh /= 0      )   CALL dom_wri      ! Create a domain file 
  • branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90

    r6943 r6952  
    261261      IF( ln_dust .AND. ln_solub ) THEN               ;  ll_solub = .TRUE. 
    262262      ELSE                                            ;  ll_solub = .FALSE. 
    263       ENDIF 
    264  
    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 (:,:) = fsdept(:,:,1) 
    270263      ENDIF 
    271264 
  • branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/TOP_SRC/TRP/trcadv.F90

    r5385 r6952  
    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         IF( lk_traldf_eiv .AND. .NOT. ln_traldf_grif )   &  ! add the eiv transport (if necessary) 
     115            &              CALL tra_adv_eiv( kt, nittrc000, zun, zvn, zwn, 'TRC' ) 
     116         ! 
     117         IF( ln_mle    )   CALL tra_adv_mle( kt, nittrc000, zun, zvn, zwn, 'TRC' )    ! add the mle transport (if necessary) 
     118         ! 
    102119      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) 
    112120      ! 
    113121      SELECT CASE ( nadv )                            !==  compute advection trend and add it to general trend  ==! 
  • branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/TOP_SRC/TRP/trcnxt.F90

    r6204 r6952  
    4545   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::   r2dt 
    4646 
     47   !! * Substitutions 
     48#  include "domzgr_substitute.h90" 
     49#  include "vectopt_loop_substitute.h90" 
    4750   !!---------------------------------------------------------------------- 
    4851   !! NEMO/TOP 3.3 , NEMO Consortium (2010) 
     
    5053   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    5154   !!---------------------------------------------------------------------- 
     55 
    5256CONTAINS 
    5357 
     
    136140         !                                               
    137141      ELSE 
    138          ! Leap-Frog + Asselin filter time stepping 
    139          IF( lk_vvl ) THEN   ;   CALL tra_nxt_vvl( kt, nittrc000, rdttrc, 'TRC', trb, trn, tra,      & 
    140            &                                                                sbc_trc, sbc_trc_b, jptra )      ! variable volume level (vvl)  
    141          ELSE                ;   CALL tra_nxt_fix( kt, nittrc000,         'TRC', trb, trn, tra, jptra )      ! fixed    volume level  
     142         IF( .NOT. lk_offline ) THEN ! Leap-Frog + Asselin filter time stepping 
     143            IF( lk_vvl ) THEN   ;   CALL tra_nxt_vvl( kt, nittrc000, rdttrc, 'TRC', trb, trn, tra,      & 
     144              &                                                                sbc_trc, sbc_trc_b, jptra )      ! variable volume level (vvl)  
     145            ELSE                ;   CALL tra_nxt_fix( kt, nittrc000,         'TRC', trb, trn, tra, jptra )      ! fixed    volume level  
     146            ENDIF 
     147         ELSE 
     148                                   CALL trc_nxt_off( kt )       ! offline  
    142149         ENDIF 
    143150      ENDIF 
     
    165172   END SUBROUTINE trc_nxt 
    166173 
     174   SUBROUTINE trc_nxt_off( kt ) 
     175      !!---------------------------------------------------------------------- 
     176      !!                   ***  ROUTINE tra_nxt_vvl  *** 
     177      !! 
     178      !! ** Purpose :   Time varying volume: apply the Asselin time filter   
     179      !!                and swap the tracer fields. 
     180      !!  
     181      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields. 
     182      !!              - save in (ta,sa) a thickness weighted average over the three  
     183      !!             time levels which will be used to compute rdn and thus the semi- 
     184      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T) 
     185      !!              - swap tracer fields to prepare the next time_step. 
     186      !!                This can be summurized for tempearture as: 
     187      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T 
     188      !!                  /( e3t_n    + rbcp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )    
     189      !!             ztm = 0                                                       otherwise 
     190      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) 
     191      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] ) 
     192      !!             tn  = ta  
     193      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call) 
     194      !! 
     195      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step 
     196      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T) 
     197      !!---------------------------------------------------------------------- 
     198      INTEGER         , INTENT(in   )                               ::  kt       ! ocean time-step index 
     199      !!      
     200      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices 
     201      REAL(wp) ::   zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar 
     202      REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      - 
     203      !!---------------------------------------------------------------------- 
     204      ! 
     205      IF( kt == nittrc000 )  THEN 
     206         IF(lwp) WRITE(numout,*) 
     207         IF(lwp) WRITE(numout,*) 'trc_nxt_off : time stepping' 
     208         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
     209      ENDIF 
     210      ! 
     211      IF(lwp) WRITE(numout,*) 'trc_nxt_off : time stepping kt = ', kt 
     212      IF(lwp) WRITE(numout,*)  
     213 
     214      DO jn = 1, jptra       
     215         DO jk = 1, jpkm1 
     216            zfact1 = atfp * rdttrc(jk) 
     217            zfact2 = zfact1 / rau0 
     218            DO jj = 1, jpj 
     219               DO ji = 1, jpi 
     220                  ze3t_b = fse3t_b(ji,jj,jk) 
     221                  ze3t_n = fse3t_n(ji,jj,jk) 
     222                  ze3t_a = fse3t_a(ji,jj,jk) 
     223                  !                                         ! tracer content at Before, now and after 
     224                  ztc_b  = trb(ji,jj,jk,jn) * ze3t_b 
     225                  ztc_n  = trn(ji,jj,jk,jn) * ze3t_n 
     226                  ztc_a  = tra(ji,jj,jk,jn) * ze3t_a 
     227                  ! 
     228                  ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b 
     229                  ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b 
     230                  ! 
     231                  ze3t_f = ze3t_n + atfp * ze3t_d 
     232                  ztc_f  = ztc_n  + atfp * ztc_d 
     233                  ! 
     234                  IF( lk_vvl .AND. jk == mikt(ji,jj) ) THEN           ! first level  
     235                     ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj)      - emp(ji,jj)   )  
     236                     ztc_f  = ztc_f  - zfact1 * ( sbc_trc(ji,jj,jn) - sbc_trc_b(ji,jj,jn) ) 
     237                  ENDIF 
     238 
     239                  ze3t_f = 1.e0 / ze3t_f 
     240                  trb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered 
     241                  trn(ji,jj,jk,jn) = tra(ji,jj,jk,jn)     ! ptn <-- pta 
     242                  ! 
     243               END DO 
     244            END DO 
     245         END DO 
     246         !  
     247      END DO 
     248      ! 
     249   END SUBROUTINE trc_nxt_off 
    167250#else 
    168251   !!---------------------------------------------------------------------- 
  • branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/TOP_SRC/TRP/trcrad.F90

    r4990 r6952  
    6464      IF( lk_c14b    )   CALL trc_rad_sms( kt, trb, trn, jp_c14b0, jp_c14b1              )  ! bomb C14 
    6565      IF( lk_pisces  )   CALL trc_rad_sms( kt, trb, trn, jp_pcs0 , jp_pcs1, cpreserv='Y' )  ! PISCES model 
    66       IF( lk_my_trc  )   CALL trc_rad_sms( kt, trb, trn, jp_myt0 , jp_myt1              )  ! MY_TRC model 
     66      IF( lk_my_trc  )   CALL trc_rad_sms( kt, trb, trn, jp_myt0 , jp_myt1, cpreserv='Y' )  ! MY_TRC model 
    6767 
    6868      ! 
  • branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90

    r6941 r6952  
    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  ! online coupling with vvl 
    131131         zsfx(:,:) = 0._wp 
    132132      ELSE                                      ! online coupling free surface or offline with free surface 
  • branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/TOP_SRC/trcnam.F90

    r6204 r6952  
    331331      !!--------------------------------------------------------------------- 
    332332 
    333       IF(lwp) WRITE(numout,*)  
    334       IF(lwp) WRITE(numout,*) 'trc_nam_dia : read the passive tracer diagnostics options' 
    335       IF(lwp) WRITE(numout,*) '~~~~~~~' 
    336  
    337333      IF(lwp) WRITE(numout,*) 
    338334      IF(lwp) WRITE(numout,*) 'trc_nam_dia : read the passive tracer diagnostics options' 
  • branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/TOP_SRC/trcstp.F90

    r6941 r6952  
    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.