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 13463 for NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OFF/dtadyn.F90 – NEMO

Ignore:
Timestamp:
2020-09-14T17:40:34+02:00 (4 years ago)
Author:
andmirek
Message:

Ticket #2195:update to trunk 13461

Location:
NEMO/branches/2019/dev_r11351_fldread_with_XIOS
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11351_fldread_with_XIOS

    • Property svn:externals
      •  

        old new  
        33^/utils/build/mk@HEAD         mk 
        44^/utils/tools@HEAD            tools 
        5 ^/vendors/AGRIF/dev@HEAD      ext/AGRIF 
         5^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 
         8 
         9# SETTE 
         10^/utils/CI/sette@13382        sette 
  • NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OFF/dtadyn.F90

    r10425 r13463  
    1313   !!             3.3  ! 2010-11 (C. Ethe) Full reorganization of the off-line: phasing with the on-line 
    1414   !!             3.4  ! 2011-05 (C. Ethe) Use of fldread 
     15   !!             4.1  ! 2019-08 (A. Coward, D. Storkey) split dta_dyn_sf_swp -> dta_dyn_sf_atf and dta_dyn_sf_interp 
    1516   !!---------------------------------------------------------------------- 
    1617 
     
    2223   USE c1d             ! 1D configuration: lk_c1d 
    2324   USE dom_oce         ! ocean domain: variables 
     25#if ! defined key_qco  
    2426   USE domvvl          ! variable volume 
     27#else 
     28   USE domqco 
     29#endif 
    2530   USE zdf_oce         ! ocean vertical physics: variables 
    2631   USE sbc_oce         ! surface module: variables 
     
    4651   PRIVATE 
    4752 
    48    PUBLIC   dta_dyn_init       ! called by opa.F90 
    49    PUBLIC   dta_dyn            ! called by step.F90 
    50    PUBLIC   dta_dyn_sed_init   ! called by opa.F90 
    51    PUBLIC   dta_dyn_sed        ! called by step.F90 
    52    PUBLIC   dta_dyn_swp        ! called by step.F90 
     53   PUBLIC   dta_dyn_init       ! called by nemo_init 
     54   PUBLIC   dta_dyn            ! called by nemo_gcm 
     55   PUBLIC   dta_dyn_sed_init   ! called by nemo_init 
     56   PUBLIC   dta_dyn_sed        ! called by nemo_gcm 
     57   PUBLIC   dta_dyn_atf        ! called by nemo_gcm 
     58#if ! defined key_qco 
     59   PUBLIC   dta_dyn_sf_interp  ! called by nemo_gcm 
     60#endif 
    5361 
    5462   CHARACTER(len=100) ::   cn_dir          !: Root directory for location of ssr files 
     
    6371   INTEGER  , SAVE      ::   jf_uwd         ! index of u-transport 
    6472   INTEGER  , SAVE      ::   jf_vwd         ! index of v-transport 
    65    INTEGER  , SAVE      ::   jf_wwd         ! index of v-transport 
     73   INTEGER  , SAVE      ::   jf_wwd         ! index of w-transport 
    6674   INTEGER  , SAVE      ::   jf_avt         ! index of Kz 
    6775   INTEGER  , SAVE      ::   jf_mld         ! index of mixed layer deptht 
     
    8795   INTEGER, SAVE  :: nprevrec, nsecdyn 
    8896 
     97   !! * Substitutions 
     98#  include "do_loop_substitute.h90" 
    8999   !!---------------------------------------------------------------------- 
    90100   !! NEMO/OFF 4.0 , NEMO Consortium (2018) 
     
    94104CONTAINS 
    95105 
    96    SUBROUTINE dta_dyn( kt ) 
     106   SUBROUTINE dta_dyn( kt, Kbb, Kmm, Kaa ) 
    97107      !!---------------------------------------------------------------------- 
    98108      !!                  ***  ROUTINE dta_dyn  *** 
     
    105115      !!             - interpolates data if needed 
    106116      !!---------------------------------------------------------------------- 
    107       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     117      INTEGER, INTENT(in) ::   kt             ! ocean time-step index 
     118      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa  ! ocean time level indices 
    108119      ! 
    109120      INTEGER             ::   ji, jj, jk 
     
    117128      ! 
    118129      IF( kt == nit000 ) THEN    ;    nprevrec = 0 
    119       ELSE                       ;    nprevrec = sf_dyn(jf_tem)%nrec_a(2) 
     130      ELSE                       ;    nprevrec = sf_dyn(jf_tem)%nrec(2,sf_dyn(jf_tem)%naa) 
    120131      ENDIF 
    121132      CALL fld_read( kt, 1, sf_dyn )      !=  read data at kt time step   ==! 
    122133      ! 
    123       IF( l_ldfslp .AND. .NOT.lk_c1d )   CALL  dta_dyn_slp( kt )    ! Computation of slopes 
    124       ! 
    125       tsn(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:)  * tmask(:,:,:)    ! temperature 
    126       tsn(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:)  * tmask(:,:,:)    ! salinity 
     134      IF( l_ldfslp .AND. .NOT.lk_c1d )   CALL  dta_dyn_slp( kt, Kbb, Kmm )    ! Computation of slopes 
     135      ! 
     136      ts(:,:,:,jp_tem,Kmm) = sf_dyn(jf_tem)%fnow(:,:,:)  * tmask(:,:,:)    ! temperature 
     137      ts(:,:,:,jp_sal,Kmm) = sf_dyn(jf_sal)%fnow(:,:,:)  * tmask(:,:,:)    ! salinity 
    127138      wndm(:,:)         = sf_dyn(jf_wnd)%fnow(:,:,1)  * tmask(:,:,1)    ! wind speed - needed for gas exchange 
    128139      fmmflx(:,:)       = sf_dyn(jf_fmf)%fnow(:,:,1)  * tmask(:,:,1)    ! downward salt flux (v3.5+) 
     
    132143      IF( ln_dynrnf ) THEN  
    133144         rnf (:,:)      = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask(:,:,1)    ! E-P 
    134          IF( ln_dynrnf_depth .AND. .NOT. ln_linssh )    CALL  dta_dyn_hrnf 
    135       ENDIF 
    136       ! 
    137       un(:,:,:)        = sf_dyn(jf_uwd)%fnow(:,:,:) * umask(:,:,:)    ! effective u-transport 
    138       vn(:,:,:)        = sf_dyn(jf_vwd)%fnow(:,:,:) * vmask(:,:,:)    ! effective v-transport 
    139       wn(:,:,:)        = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:)    ! effective v-transport 
     145         IF( ln_dynrnf_depth .AND. .NOT. ln_linssh )    CALL  dta_dyn_hrnf(Kmm) 
     146      ENDIF 
     147      ! 
     148      uu(:,:,:,Kmm)        = sf_dyn(jf_uwd)%fnow(:,:,:) * umask(:,:,:)    ! effective u-transport 
     149      vv(:,:,:,Kmm)        = sf_dyn(jf_vwd)%fnow(:,:,:) * vmask(:,:,:)    ! effective v-transport 
     150      ww(:,:,:)        = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:)    ! effective v-transport 
    140151      ! 
    141152      IF( .NOT.ln_linssh ) THEN 
     
    144155         emp_b  (:,:)   = sf_dyn(jf_empb)%fnow(:,:,1) * tmask(:,:,1)    ! E-P 
    145156         zemp   (:,:)   = ( 0.5_wp * ( emp(:,:) + emp_b(:,:) ) + rnf(:,:) + fwbcorr ) * tmask(:,:,1) 
    146          CALL dta_dyn_ssh( kt, zhdivtr, sshb, zemp, ssha, e3t_a(:,:,:) )  !=  ssh, vertical scale factor & vertical transport 
     157#if defined key_qco 
     158         CALL dta_dyn_ssh( kt, zhdivtr, ssh(:,:,Kbb), zemp, ssh(:,:,Kaa) ) 
     159         CALL dom_qco_r3c( ssh(:,:,Kaa), r3t(:,:,Kaa), r3u(:,:,Kaa), r3v(:,:,Kaa) ) 
     160#else 
     161         CALL dta_dyn_ssh( kt, zhdivtr, ssh(:,:,Kbb), zemp, ssh(:,:,Kaa), e3t(:,:,:,Kaa) )  !=  ssh, vertical scale factor 
     162#endif 
    147163         DEALLOCATE( zemp , zhdivtr ) 
    148164         !                                           Write in the tracer restart file 
     
    152168            IF(lwp) WRITE(numout,*) 'dta_dyn_ssh : ssh field written in tracer restart file at it= ', kt,' date= ', ndastp 
    153169            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    154             CALL iom_rstput( kt, nitrst, numrtw, 'sshn', ssha ) 
    155             CALL iom_rstput( kt, nitrst, numrtw, 'sshb', sshn ) 
     170            CALL iom_rstput( kt, nitrst, numrtw, 'sshn', ssh(:,:,Kaa) ) 
     171            CALL iom_rstput( kt, nitrst, numrtw, 'sshb', ssh(:,:,Kmm) ) 
    156172         ENDIF 
    157173      ENDIF 
    158174      ! 
    159       CALL eos    ( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 
    160       CALL eos_rab( tsn, rab_n )       ! now    local thermal/haline expension ratio at T-points 
    161       CALL bn2    ( tsn, rab_n, rn2 ) ! before Brunt-Vaisala frequency need for zdfmxl 
    162  
    163       rn2b(:,:,:) = rn2(:,:,:)         ! need for zdfmxl 
    164       CALL zdf_mxl( kt )                                                   ! In any case, we need mxl 
     175      CALL eos    ( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 
     176      CALL eos_rab( ts(:,:,:,:,Kmm), rab_n, Kmm )       ! now    local thermal/haline expension ratio at T-points 
     177      CALL bn2    ( ts(:,:,:,:,Kmm), rab_n, rn2, Kmm ) ! before Brunt-Vaisala frequency need for zdfmxl 
     178 
     179      rn2b(:,:,:) = rn2(:,:,:)         ! needed for zdfmxl 
     180      CALL zdf_mxl( kt, Kmm )          ! In any case, we need mxl 
    165181      ! 
    166182      hmld(:,:)       = sf_dyn(jf_mld)%fnow(:,:,1) * tmask(:,:,1)    ! mixed layer depht 
     
    174190      ! 
    175191      ! 
    176       CALL eos( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 
    177       ! 
    178       IF(ln_ctl) THEN                  ! print control 
    179          CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' tn      - : ', mask1=tmask,  kdim=jpk   ) 
    180          CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_sal), clinfo1=' sn      - : ', mask1=tmask,  kdim=jpk   ) 
    181          CALL prt_ctl(tab3d_1=un               , clinfo1=' un      - : ', mask1=umask,  kdim=jpk   ) 
    182          CALL prt_ctl(tab3d_1=vn               , clinfo1=' vn      - : ', mask1=vmask,  kdim=jpk   ) 
    183          CALL prt_ctl(tab3d_1=wn               , clinfo1=' wn      - : ', mask1=tmask,  kdim=jpk   ) 
     192      CALL eos( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 
     193      ! 
     194      IF(sn_cfctl%l_prtctl) THEN                 ! print control 
     195         CALL prt_ctl(tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' tn      - : ', mask1=tmask,  kdim=jpk   ) 
     196         CALL prt_ctl(tab3d_1=ts(:,:,:,jp_sal,Kmm), clinfo1=' sn      - : ', mask1=tmask,  kdim=jpk   ) 
     197         CALL prt_ctl(tab3d_1=uu(:,:,:,Kmm)               , clinfo1=' uu(:,:,:,Kmm)      - : ', mask1=umask,  kdim=jpk   ) 
     198         CALL prt_ctl(tab3d_1=vv(:,:,:,Kmm)               , clinfo1=' vv(:,:,:,Kmm)      - : ', mask1=vmask,  kdim=jpk   ) 
     199         CALL prt_ctl(tab3d_1=ww               , clinfo1=' ww      - : ', mask1=tmask,  kdim=jpk   ) 
    184200         CALL prt_ctl(tab3d_1=avt              , clinfo1=' kz      - : ', mask1=tmask,  kdim=jpk   ) 
    185201         CALL prt_ctl(tab3d_1=uslp             , clinfo1=' slp  - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk) 
     
    192208 
    193209 
    194    SUBROUTINE dta_dyn_init 
     210   SUBROUTINE dta_dyn_init( Kbb, Kmm, Kaa ) 
    195211      !!---------------------------------------------------------------------- 
    196212      !!                  ***  ROUTINE dta_dyn_init  *** 
     
    199215      !! ** Method  : - read the data namdta_dyn namelist 
    200216      !!---------------------------------------------------------------------- 
     217      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa  ! ocean time level indices 
     218      ! 
    201219      INTEGER  :: ierr, ierr0, ierr1, ierr2, ierr3   ! return error code 
    202220      INTEGER  :: ifpr                               ! dummy loop indice 
     
    225243      !!---------------------------------------------------------------------- 
    226244      ! 
    227       REWIND( numnam_ref )              ! Namelist namdta_dyn in reference namelist : Offline: init. of dynamical data 
    228245      READ  ( numnam_ref, namdta_dyn, IOSTAT = ios, ERR = 901) 
    229 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdta_dyn in reference namelist', lwp ) 
    230       REWIND( numnam_cfg )              ! Namelist namdta_dyn in configuration namelist : Offline: init. of dynamical data 
     246901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdta_dyn in reference namelist' ) 
    231247      READ  ( numnam_cfg, namdta_dyn, IOSTAT = ios, ERR = 902 ) 
    232 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdta_dyn in configuration namelist', lwp ) 
     248902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdta_dyn in configuration namelist' ) 
    233249      IF(lwm) WRITE ( numond, namdta_dyn ) 
    234250      !                                         ! store namelist information in an array 
     
    278294      !                                         ! fill sf with slf_i and control print 
    279295      CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' ) 
     296      sf_dyn(jf_uwd)%cltype = 'U'   ;   sf_dyn(jf_uwd)%zsgn = -1._wp   
     297      sf_dyn(jf_vwd)%cltype = 'V'   ;   sf_dyn(jf_vwd)%zsgn = -1._wp   
     298      ! 
     299      IF( ln_trabbl ) THEN 
     300         sf_dyn(jf_ubl)%cltype = 'U'   ;   sf_dyn(jf_ubl)%zsgn =  1._wp   
     301         sf_dyn(jf_vbl)%cltype = 'V'   ;   sf_dyn(jf_vbl)%zsgn =  1._wp   
     302      END IF 
    280303      ! 
    281304      ! Open file for each variable to get his number of dimension 
    282305      DO ifpr = 1, jfld 
    283          CALL fld_clopn( sf_dyn(ifpr), nyear, nmonth, nday ) 
     306         CALL fld_def( sf_dyn(ifpr) ) 
     307         CALL iom_open( sf_dyn(ifpr)%clname, sf_dyn(ifpr)%num ) 
    284308         idv   = iom_varid( sf_dyn(ifpr)%num , slf_d(ifpr)%clvar )        ! id of the variable sdjf%clvar 
    285309         idimv = iom_file ( sf_dyn(ifpr)%num )%ndims(idv)                 ! number of dimension for variable sdjf%clvar 
    286          IF( sf_dyn(ifpr)%num /= 0 )   CALL iom_close( sf_dyn(ifpr)%num ) ! close file if already open 
     310         CALL iom_close( sf_dyn(ifpr)%num )                              ! close file if already open 
    287311         ierr1=0 
    288312         IF( idimv == 3 ) THEN    ! 2D variable 
     
    312336        IF( .NOT. sf_dyn(jf_uwd)%ln_clim .AND. ln_rsttr .AND.    &                     ! Restart: read in restart file 
    313337           iom_varid( numrtr, 'sshn', ldstop = .FALSE. ) > 0 ) THEN 
    314            IF(lwp) WRITE(numout,*) ' sshn forcing fields read in the restart file for initialisation' 
    315            CALL iom_get( numrtr, jpdom_autoglo, 'sshn', sshn(:,:)   ) 
    316            CALL iom_get( numrtr, jpdom_autoglo, 'sshb', sshb(:,:)   ) 
     338           IF(lwp) WRITE(numout,*) ' ssh(:,:,Kmm) forcing fields read in the restart file for initialisation' 
     339           CALL iom_get( numrtr, jpdom_auto, 'sshn', ssh(:,:,Kmm)   ) 
     340           CALL iom_get( numrtr, jpdom_auto, 'sshb', ssh(:,:,Kbb)   ) 
    317341        ELSE 
    318            IF(lwp) WRITE(numout,*) ' sshn forcing fields read in the restart file for initialisation' 
     342           IF(lwp) WRITE(numout,*) ' ssh(:,:,Kmm) forcing fields read in the restart file for initialisation' 
    319343           CALL iom_open( 'restart', inum ) 
    320            CALL iom_get( inum, jpdom_autoglo, 'sshn', sshn(:,:)   ) 
    321            CALL iom_get( inum, jpdom_autoglo, 'sshb', sshb(:,:)   ) 
     344           CALL iom_get( inum, jpdom_auto, 'sshn', ssh(:,:,Kmm)   ) 
     345           CALL iom_get( inum, jpdom_auto, 'sshb', ssh(:,:,Kbb)   ) 
    322346           CALL iom_close( inum )                                        ! close file 
    323347        ENDIF 
    324348        ! 
     349#if defined key_qco 
     350        CALL dom_qco_r3c( ssh(:,:,Kbb), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb) ) 
     351        CALL dom_qco_r3c( ssh(:,:,Kmm), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm) ) 
     352#else 
    325353        DO jk = 1, jpkm1 
    326            e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) ) 
     354           e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( 1._wp + ssh(:,:,Kmm) * r1_ht_0(:,:) * tmask(:,:,jk) ) 
    327355        ENDDO 
    328         e3t_a(:,:,jpk) = e3t_0(:,:,jpk) 
     356        e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk) 
    329357 
    330358        ! Horizontal scale factor interpolations 
    331359        ! -------------------------------------- 
    332         CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 
    333         CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 
     360        CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 
     361        CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 
    334362 
    335363        ! Vertical scale factor interpolations 
    336364        ! ------------------------------------ 
    337         CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n(:,:,:), 'W' ) 
    338    
    339         e3t_b(:,:,:)  = e3t_n(:,:,:) 
    340         e3u_b(:,:,:)  = e3u_n(:,:,:) 
    341         e3v_b(:,:,:)  = e3v_n(:,:,:) 
     365        CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w(:,:,:,Kmm), 'W' ) 
     366!!gm this should be computed from ssh(Kbb)   
     367        e3t(:,:,:,Kbb)  = e3t(:,:,:,Kmm) 
     368        e3u(:,:,:,Kbb)  = e3u(:,:,:,Kmm) 
     369        e3v(:,:,:,Kbb)  = e3v(:,:,:,Kmm) 
    342370 
    343371        ! t- and w- points depth 
    344372        ! ---------------------- 
    345         gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 
    346         gdepw_n(:,:,1) = 0.0_wp 
    347  
    348         DO jk = 2, jpk 
    349            DO jj = 1,jpj 
    350               DO ji = 1,jpi 
    351                 !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere 
    352                 !    tmask = wmask, ie everywhere expect at jk = mikt 
    353                                                                    ! 1 for jk = 
    354                                                                    ! mikt 
    355                  zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
    356                  gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 
    357                  gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  & 
    358                      &                + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)) 
    359               END DO 
    360            END DO 
    361         END DO 
    362  
    363         gdept_b(:,:,:) = gdept_n(:,:,:) 
    364         gdepw_b(:,:,:) = gdepw_n(:,:,:) 
     373        gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 
     374        gdepw(:,:,1,Kmm) = 0.0_wp 
     375 
     376        DO_3D( 1, 1, 1, 1, 2, jpk ) 
     377          !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere 
     378          !    tmask = wmask, ie everywhere expect at jk = mikt 
     379                                                             ! 1 for jk = 
     380                                                             ! mikt 
     381           zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
     382           gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
     383           gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  & 
     384               &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm)) 
     385        END_3D 
     386 
     387        gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm) 
     388        gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm) 
    365389        ! 
    366390      ENDIF 
     391#endif 
    367392      ! 
    368393      IF( ln_dynrnf .AND. ln_dynrnf_depth ) THEN       ! read depht over which runoffs are distributed 
     
    370395         IF(lwp) WRITE(numout,*) ' read in the file depht over which runoffs are distributed' 
    371396         CALL iom_open ( "runoffs", inum )                           ! open file 
    372          CALL iom_get  ( inum, jpdom_data, 'rodepth', h_rnf )   ! read the river mouth array 
     397         CALL iom_get  ( inum, jpdom_global, 'rodepth', h_rnf )   ! read the river mouth array 
    373398         CALL iom_close( inum )                                        ! close file 
    374399         ! 
    375400         nk_rnf(:,:) = 0                               ! set the number of level over which river runoffs are applied 
    376          DO jj = 1, jpj 
    377             DO ji = 1, jpi 
    378                IF( h_rnf(ji,jj) > 0._wp ) THEN 
    379                   jk = 2 
    380                   DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1 
    381                   END DO 
    382                   nk_rnf(ji,jj) = jk 
    383                ELSEIF( h_rnf(ji,jj) == -1._wp   ) THEN   ;  nk_rnf(ji,jj) = 1 
    384                ELSEIF( h_rnf(ji,jj) == -999._wp ) THEN   ;  nk_rnf(ji,jj) = mbkt(ji,jj) 
    385                ELSE 
    386                   CALL ctl_stop( 'sbc_rnf_init: runoff depth not positive, and not -999 or -1, rnf value in file fort.999'  ) 
    387                   WRITE(999,*) 'ji, jj, h_rnf(ji,jj) :', ji, jj, h_rnf(ji,jj) 
    388                ENDIF 
     401         DO_2D( 1, 1, 1, 1 ) 
     402            IF( h_rnf(ji,jj) > 0._wp ) THEN 
     403               jk = 2 
     404               DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1 
     405               END DO 
     406               nk_rnf(ji,jj) = jk 
     407            ELSEIF( h_rnf(ji,jj) == -1._wp   ) THEN   ;  nk_rnf(ji,jj) = 1 
     408            ELSEIF( h_rnf(ji,jj) == -999._wp ) THEN   ;  nk_rnf(ji,jj) = mbkt(ji,jj) 
     409            ELSE 
     410               CALL ctl_stop( 'sbc_rnf_init: runoff depth not positive, and not -999 or -1, rnf value in file fort.999'  ) 
     411               WRITE(999,*) 'ji, jj, h_rnf(ji,jj) :', ji, jj, h_rnf(ji,jj) 
     412            ENDIF 
     413         END_2D 
     414!!st pourquoi on n'utilise pas le gde3w ici plutôt que de faire une boucle ?  
     415         DO_2D( 1, 1, 1, 1 ) 
     416            h_rnf(ji,jj) = 0._wp 
     417            DO jk = 1, nk_rnf(ji,jj) 
     418               h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm) 
    389419            END DO 
    390          END DO 
    391          DO jj = 1, jpj                                ! set the associated depth 
    392             DO ji = 1, jpi 
    393                h_rnf(ji,jj) = 0._wp 
    394                DO jk = 1, nk_rnf(ji,jj) 
    395                   h_rnf(ji,jj) = h_rnf(ji,jj) + e3t_n(ji,jj,jk) 
    396                END DO 
    397             END DO 
    398          END DO 
     420         END_2D 
    399421      ELSE                                       ! runoffs applied at the surface 
    400422         nk_rnf(:,:) = 1 
    401          h_rnf (:,:) = e3t_n(:,:,1) 
     423         h_rnf (:,:) = e3t(:,:,1,Kmm) 
    402424      ENDIF 
    403425      nkrnf_max = MAXVAL( nk_rnf(:,:) ) 
     
    411433      IF(lwp) WRITE(numout,*) ' ' 
    412434      ! 
    413       CALL dta_dyn( nit000 ) 
     435      CALL dta_dyn( nit000, Kbb, Kmm, Kaa ) 
    414436      ! 
    415437   END SUBROUTINE dta_dyn_init 
    416438 
    417    SUBROUTINE dta_dyn_sed( kt ) 
     439    
     440   SUBROUTINE dta_dyn_sed( kt, Kmm ) 
    418441      !!---------------------------------------------------------------------- 
    419442      !!                  ***  ROUTINE dta_dyn  *** 
     
    427450      !!---------------------------------------------------------------------- 
    428451      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     452      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index 
    429453      ! 
    430454      !!---------------------------------------------------------------------- 
     
    435459      ! 
    436460      IF( kt == nit000 ) THEN    ;    nprevrec = 0 
    437       ELSE                       ;    nprevrec = sf_dyn(jf_tem)%nrec_a(2) 
     461      ELSE                       ;    nprevrec = sf_dyn(jf_tem)%nrec(2,sf_dyn(jf_tem)%naa) 
    438462      ENDIF 
    439463      CALL fld_read( kt, 1, sf_dyn )      !=  read data at kt time step   ==! 
    440464      ! 
    441       tsn(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:)  * tmask(:,:,:)    ! temperature 
    442       tsn(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:)  * tmask(:,:,:)    ! salinity 
    443       ! 
    444       CALL eos    ( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 
    445  
    446       IF(ln_ctl) THEN                  ! print control 
    447          CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' tn      - : ', mask1=tmask,  kdim=jpk   ) 
    448          CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_sal), clinfo1=' sn      - : ', mask1=tmask,  kdim=jpk   ) 
     465      ts(:,:,:,jp_tem,Kmm) = sf_dyn(jf_tem)%fnow(:,:,:)  * tmask(:,:,:)    ! temperature 
     466      ts(:,:,:,jp_sal,Kmm) = sf_dyn(jf_sal)%fnow(:,:,:)  * tmask(:,:,:)    ! salinity 
     467      ! 
     468      CALL eos    ( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 
     469 
     470      IF(sn_cfctl%l_prtctl) THEN                     ! print control 
     471         CALL prt_ctl(tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' tn      - : ', mask1=tmask,  kdim=jpk   ) 
     472         CALL prt_ctl(tab3d_1=ts(:,:,:,jp_sal,Kmm), clinfo1=' sn      - : ', mask1=tmask,  kdim=jpk   ) 
    449473      ENDIF 
    450474      ! 
     
    454478 
    455479 
    456    SUBROUTINE dta_dyn_sed_init 
     480   SUBROUTINE dta_dyn_sed_init( Kmm ) 
    457481      !!---------------------------------------------------------------------- 
    458482      !!                  ***  ROUTINE dta_dyn_init  *** 
     
    461485      !! ** Method  : - read the data namdta_dyn namelist 
    462486      !!---------------------------------------------------------------------- 
     487      INTEGER, INTENT( in ) :: Kmm                   ! ocean time level index 
     488      ! 
    463489      INTEGER  :: ierr, ierr0, ierr1, ierr2, ierr3   ! return error code 
    464490      INTEGER  :: ifpr                               ! dummy loop indice 
     
    475501      !!---------------------------------------------------------------------- 
    476502      ! 
    477       REWIND( numnam_ref )              ! Namelist namdta_dyn in reference namelist : Offline: init. of dynamical data 
    478503      READ  ( numnam_ref, namdta_dyn, IOSTAT = ios, ERR = 901) 
    479 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdta_dyn in reference namelist', lwp ) 
    480       REWIND( numnam_cfg )              ! Namelist namdta_dyn in configuration namelist : Offline: init. of dynamical data 
     504901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdta_dyn in reference namelist' ) 
    481505      READ  ( numnam_cfg, namdta_dyn, IOSTAT = ios, ERR = 902 ) 
    482 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdta_dyn in configuration namelist', lwp ) 
     506902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdta_dyn in configuration namelist' ) 
    483507      IF(lwm) WRITE ( numond, namdta_dyn ) 
    484508      !                                         ! store namelist information in an array 
     
    508532      ! Open file for each variable to get his number of dimension 
    509533      DO ifpr = 1, jfld 
    510          CALL fld_clopn( sf_dyn(ifpr), nyear, nmonth, nday ) 
     534         CALL fld_def( sf_dyn(ifpr) ) 
     535         CALL iom_open( sf_dyn(ifpr)%clname, sf_dyn(ifpr)%num ) 
    511536         idv   = iom_varid( sf_dyn(ifpr)%num , slf_d(ifpr)%clvar )        ! id of the variable sdjf%clvar 
    512537         idimv = iom_file ( sf_dyn(ifpr)%num )%ndims(idv)                 ! number of dimension for variable sdjf%clvar 
    513          IF( sf_dyn(ifpr)%num /= 0 )   CALL iom_close( sf_dyn(ifpr)%num ) ! close file if already open 
     538         CALL iom_close( sf_dyn(ifpr)%num )                              ! close file if already open 
    514539         ierr1=0 
    515540         IF( idimv == 3 ) THEN    ! 2D variable 
     
    525550      END DO 
    526551      ! 
    527       CALL dta_dyn_sed( nit000 ) 
     552      CALL dta_dyn_sed( nit000, Kmm ) 
    528553      ! 
    529554   END SUBROUTINE dta_dyn_sed_init 
    530555 
    531    SUBROUTINE dta_dyn_swp( kt ) 
     556    
     557   SUBROUTINE dta_dyn_atf( kt, Kbb, Kmm, Kaa ) 
    532558     !!--------------------------------------------------------------------- 
    533559      !!                    ***  ROUTINE dta_dyn_swp  *** 
    534560      !! 
    535       !! ** Purpose :   Swap and the data and compute the vertical scale factor  
    536       !!              at U/V/W pointand the depht 
    537       !!--------------------------------------------------------------------- 
    538       INTEGER, INTENT(in) :: kt       ! time step 
     561      !! ** Purpose :   Asselin time filter of now SSH 
     562      !!--------------------------------------------------------------------- 
     563      INTEGER, INTENT(in) :: kt             ! time step 
     564      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa  ! ocean time level indices 
     565      ! 
     566      !!--------------------------------------------------------------------- 
     567 
     568      IF( kt == nit000 ) THEN 
     569         IF(lwp) WRITE(numout,*) 
     570         IF(lwp) WRITE(numout,*) 'dta_dyn_atf : Asselin time filter of sea surface height' 
     571         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 
     572      ENDIF 
     573 
     574      ssh(:,:,Kmm) = ssh(:,:,Kmm) + rn_atfp * ( ssh(:,:,Kbb) - 2 * ssh(:,:,Kmm) + ssh(:,:,Kaa))   
     575 
     576      !! Do we also need to time filter e3t?? 
     577      ! 
     578   END SUBROUTINE dta_dyn_atf 
     579    
     580    
     581#if ! defined key_qco     
     582   SUBROUTINE dta_dyn_sf_interp( kt, Kmm ) 
     583      !!--------------------------------------------------------------------- 
     584      !!                    ***  ROUTINE dta_dyn_sf_interp  *** 
     585      !! 
     586      !! ** Purpose :   Calculate scale factors at U/V/W points and depths 
     587      !!                given the after e3t field 
     588      !!--------------------------------------------------------------------- 
     589      INTEGER, INTENT(in) :: kt   ! time step 
     590      INTEGER, INTENT(in) :: Kmm  ! ocean time level indices 
    539591      ! 
    540592      INTEGER             :: ji, jj, jk 
     
    542594      !!--------------------------------------------------------------------- 
    543595 
    544       IF( kt == nit000 ) THEN 
    545          IF(lwp) WRITE(numout,*) 
    546          IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height' 
    547          IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
    548       ENDIF 
    549  
    550       sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:))  ! before <-- now filtered 
    551       sshn(:,:) = ssha(:,:) 
    552  
    553       e3t_n(:,:,:) = e3t_a(:,:,:) 
    554  
    555       ! Reconstruction of all vertical scale factors at now and before time steps 
    556       ! ============================================================================= 
    557  
    558596      ! Horizontal scale factor interpolations 
    559597      ! -------------------------------------- 
    560       CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 
    561       CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 
     598      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 
     599      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 
    562600 
    563601      ! Vertical scale factor interpolations 
    564602      ! ------------------------------------ 
    565       CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' ) 
    566  
    567       e3t_b(:,:,:)  = e3t_n(:,:,:) 
    568       e3u_b(:,:,:)  = e3u_n(:,:,:) 
    569       e3v_b(:,:,:)  = e3v_n(:,:,:) 
     603      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' ) 
    570604 
    571605      ! t- and w- points depth 
    572606      ! ---------------------- 
    573       gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 
    574       gdepw_n(:,:,1) = 0.0_wp 
    575       ! 
    576       DO jk = 2, jpk 
    577          DO jj = 1,jpj 
    578             DO ji = 1,jpi 
    579                zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
    580                gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 
    581                gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  & 
    582                   &                + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)) 
    583             END DO 
    584          END DO 
    585       END DO 
    586       ! 
    587       gdept_b(:,:,:) = gdept_n(:,:,:) 
    588       gdepw_b(:,:,:) = gdepw_n(:,:,:) 
    589       ! 
    590    END SUBROUTINE dta_dyn_swp 
    591     
     607      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 
     608      gdepw(:,:,1,Kmm) = 0.0_wp 
     609      ! 
     610      DO_3D( 1, 1, 1, 1, 2, jpk ) 
     611         zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
     612         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
     613         gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  & 
     614            &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm)) 
     615      END_3D 
     616      ! 
     617   END SUBROUTINE dta_dyn_sf_interp 
     618#endif 
     619 
    592620 
    593621   SUBROUTINE dta_dyn_ssh( kt, phdivtr, psshb,  pemp, pssha, pe3ta ) 
     
    595623      !!                ***  ROUTINE dta_dyn_wzv  *** 
    596624      !!                    
    597       !! ** Purpose :   compute the after ssh (ssha) and the now vertical velocity 
     625      !! ** Purpose :   compute the after ssh (ssh(:,:,Kaa)) and the now vertical velocity 
    598626      !! 
    599627      !! ** Method  : Using the incompressibility hypothesis,  
     
    608636      !!          The boundary conditions are w=0 at the bottom (no flux) 
    609637      !! 
    610       !! ** action  :   ssha / e3t_a / wn 
     638      !! ** action  :   ssh(:,:,Kaa) / e3t(:,:,k,Kaa) / ww 
    611639      !! 
    612640      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
     
    624652      !!---------------------------------------------------------------------- 
    625653      ! 
    626       z2dt = 2._wp * rdt 
     654      z2dt = 2._wp * rn_Dt 
    627655      ! 
    628656      zhdiv(:,:) = 0._wp 
     
    631659      END DO 
    632660      !                                                ! Sea surface  elevation time-stepping 
    633       pssha(:,:) = ( psshb(:,:) - z2dt * ( r1_rau0 * pemp(:,:)  + zhdiv(:,:) ) ) * ssmask(:,:) 
    634       !                                                 !  
    635       !                                                 ! After acale factors at t-points ( z_star coordinate ) 
     661      pssha(:,:) = ( psshb(:,:) - z2dt * ( r1_rho0 * pemp(:,:)  + zhdiv(:,:) ) ) * ssmask(:,:) 
     662      ! 
     663      IF( PRESENT( pe3ta ) ) THEN                      ! After acale factors at t-points ( z_star coordinate ) 
    636664      DO jk = 1, jpkm1 
    637         pe3ta(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + pssha(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) ) 
     665            pe3ta(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + pssha(:,:) * r1_ht_0(:,:) * tmask(:,:,jk) ) 
    638666      END DO 
     667      ENDIF 
    639668      ! 
    640669   END SUBROUTINE dta_dyn_ssh 
    641670 
    642671 
    643    SUBROUTINE dta_dyn_hrnf 
     672   SUBROUTINE dta_dyn_hrnf( Kmm ) 
    644673      !!---------------------------------------------------------------------- 
    645674      !!                  ***  ROUTINE sbc_rnf  *** 
     
    654683      !!---------------------------------------------------------------------- 
    655684      !! 
    656       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    657       !!---------------------------------------------------------------------- 
    658       ! 
    659       DO jj = 1, jpj                   ! update the depth over which runoffs are distributed 
    660          DO ji = 1, jpi 
    661             h_rnf(ji,jj) = 0._wp 
    662             DO jk = 1, nk_rnf(ji,jj)                           ! recalculates h_rnf to be the depth in metres 
    663                 h_rnf(ji,jj) = h_rnf(ji,jj) + e3t_n(ji,jj,jk)   ! to the bottom of the relevant grid box 
    664             END DO 
    665         END DO 
    666       END DO 
     685      INTEGER, INTENT(in) :: Kmm  ! ocean time level index 
     686      ! 
     687      INTEGER  ::   ji, jj, jk    ! dummy loop indices 
     688      !!---------------------------------------------------------------------- 
     689      ! 
     690!!st code dupliqué même remarque que plus haut pourquoi ne pas utiliser gdepw ? 
     691      DO_2D( 1, 1, 1, 1 ) 
     692         h_rnf(ji,jj) = 0._wp 
     693         DO jk = 1, nk_rnf(ji,jj)                           ! recalculates h_rnf to be the depth in metres 
     694             h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm)   ! to the bottom of the relevant grid box 
     695         END DO 
     696      END_2D 
    667697      ! 
    668698   END SUBROUTINE dta_dyn_hrnf 
     
    670700 
    671701 
    672    SUBROUTINE dta_dyn_slp( kt ) 
     702   SUBROUTINE dta_dyn_slp( kt, Kbb, Kmm ) 
    673703      !!--------------------------------------------------------------------- 
    674704      !!                    ***  ROUTINE dta_dyn_slp  *** 
     
    678708      !!--------------------------------------------------------------------- 
    679709      INTEGER,  INTENT(in) :: kt       ! time step 
     710      INTEGER,  INTENT(in) :: Kbb, Kmm ! ocean time level indices 
    680711      ! 
    681712      INTEGER  ::   ji, jj     ! dummy loop indices 
     
    687718      !!--------------------------------------------------------------------- 
    688719      ! 
    689       IF( sf_dyn(jf_tem)%ln_tint ) THEN    ! Computes slopes (here avt is used as workspace)                        
     720      IF( sf_dyn(jf_tem)%ln_tint ) THEN    ! Computes slopes (here avt is used as workspace) 
     721         ! 
    690722         IF( kt == nit000 ) THEN 
    691723            IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt 
    692             zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,1) * tmask(:,:,:)   ! temperature 
    693             zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:)   ! salinity  
    694             avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,1) * tmask(:,:,:)   ! vertical diffusive coef. 
    695             CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
     724            zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,sf_dyn(jf_tem)%nbb) * tmask(:,:,:)   ! temperature 
     725            zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,sf_dyn(jf_sal)%nbb) * tmask(:,:,:)   ! salinity  
     726            avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,sf_dyn(jf_avt)%nbb) * tmask(:,:,:)   ! vertical diffusive coef. 
     727            CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm ) 
    696728            uslpdta (:,:,:,1) = zuslp (:,:,:)  
    697729            vslpdta (:,:,:,1) = zvslp (:,:,:)  
     
    699731            wslpjdta(:,:,:,1) = zwslpj(:,:,:)  
    700732            ! 
    701             zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:)   ! temperature 
    702             zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity  
    703             avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef. 
    704             CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
     733            zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,sf_dyn(jf_tem)%naa) * tmask(:,:,:)   ! temperature 
     734            zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,sf_dyn(jf_sal)%naa) * tmask(:,:,:)   ! salinity  
     735            avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,sf_dyn(jf_avt)%naa) * tmask(:,:,:)   ! vertical diffusive coef. 
     736            CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm ) 
    705737            uslpdta (:,:,:,2) = zuslp (:,:,:)  
    706738            vslpdta (:,:,:,2) = zvslp (:,:,:)  
     
    710742           !  
    711743           iswap = 0 
    712            IF( sf_dyn(jf_tem)%nrec_a(2) - nprevrec /= 0 )  iswap = 1 
    713            IF( nsecdyn > sf_dyn(jf_tem)%nrec_b(2) .AND. iswap == 1 )  THEN    ! read/update the after data 
     744           IF( sf_dyn(jf_tem)%nrec(2,sf_dyn(jf_tem)%naa) - nprevrec /= 0 )  iswap = 1 
     745           IF( nsecdyn > sf_dyn(jf_tem)%nrec(2,sf_dyn(jf_tem)%nbb) .AND. iswap == 1 )  THEN    ! read/update the after data 
    714746              IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt 
    715747              uslpdta (:,:,:,1) =  uslpdta (:,:,:,2)         ! swap the data 
     
    718750              wslpjdta(:,:,:,1) =  wslpjdta(:,:,:,2)  
    719751              ! 
    720               zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:)   ! temperature 
    721               zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity  
    722               avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef. 
    723               CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
     752              zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,sf_dyn(jf_tem)%naa) * tmask(:,:,:)   ! temperature 
     753              zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,sf_dyn(jf_sal)%naa) * tmask(:,:,:)   ! salinity  
     754              avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,sf_dyn(jf_avt)%naa) * tmask(:,:,:)   ! vertical diffusive coef. 
     755              CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm ) 
    724756              ! 
    725757              uslpdta (:,:,:,2) = zuslp (:,:,:)  
     
    732764      ! 
    733765      IF( sf_dyn(jf_tem)%ln_tint )  THEN 
    734          ztinta =  REAL( nsecdyn - sf_dyn(jf_tem)%nrec_b(2), wp )  & 
    735             &    / REAL( sf_dyn(jf_tem)%nrec_a(2) - sf_dyn(jf_tem)%nrec_b(2), wp ) 
     766         ztinta =  REAL( nsecdyn - sf_dyn(jf_tem)%nrec(2,sf_dyn(jf_tem)%nbb), wp )  & 
     767            &    / REAL( sf_dyn(jf_tem)%nrec(2,sf_dyn(jf_tem)%naa) - sf_dyn(jf_tem)%nrec(2,sf_dyn(jf_tem)%nbb), wp ) 
    736768         ztintb =  1. - ztinta 
    737769         IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace) 
     
    745777         zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:)   ! salinity  
    746778         avt(:,:,:)        = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:)   ! vertical diffusive coef. 
    747          CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
     779         CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm ) 
    748780         ! 
    749781         IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace) 
     
    758790 
    759791 
    760    SUBROUTINE compute_slopes( kt, pts, puslp, pvslp, pwslpi, pwslpj ) 
     792   SUBROUTINE compute_slopes( kt, pts, puslp, pvslp, pwslpi, pwslpj, Kbb, Kmm ) 
    761793      !!--------------------------------------------------------------------- 
    762794      !!                    ***  ROUTINE dta_dyn_slp  *** 
     
    770802      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pwslpi   ! zonal diapycnal slopes 
    771803      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pwslpj   ! meridional diapycnal slopes 
     804      INTEGER ,                              INTENT(in ) :: Kbb, Kmm ! ocean time level indices 
    772805      !!--------------------------------------------------------------------- 
    773806      ! 
    774807      IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace) 
    775808         CALL eos    ( pts, rhd, rhop, gdept_0(:,:,:) ) 
    776          CALL eos_rab( pts, rab_n )       ! now local thermal/haline expension ratio at T-points 
    777          CALL bn2    ( pts, rab_n, rn2  ) ! now    Brunt-Vaisala 
     809         CALL eos_rab( pts, rab_n, Kmm )       ! now local thermal/haline expension ratio at T-points 
     810         CALL bn2    ( pts, rab_n, rn2, Kmm  ) ! now    Brunt-Vaisala 
    778811 
    779812      ! Partial steps: before Horizontal DErivative 
    780813      IF( ln_zps  .AND. .NOT. ln_isfcav)                            & 
    781          &            CALL zps_hde    ( kt, jpts, pts, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient 
     814         &            CALL zps_hde    ( kt, Kmm, jpts, pts, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient 
    782815         &                                        rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
    783816      IF( ln_zps .AND.        ln_isfcav)                            & 
    784          &            CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv, gtui, gtvi, &  ! Partial steps for top cell (ISF) 
     817         &            CALL zps_hde_isf( kt, Kmm, jpts, pts, gtsu, gtsv, gtui, gtvi, &  ! Partial steps for top cell (ISF) 
    785818         &                                        rhd, gru , grv , grui, grvi )  ! of t, s, rd at the first ocean level 
    786819 
    787          rn2b(:,:,:) = rn2(:,:,:)         ! need for zdfmxl 
    788          CALL zdf_mxl( kt )            ! mixed layer depth 
    789          CALL ldf_slp( kt, rhd, rn2 )  ! slopes 
     820         rn2b(:,:,:) = rn2(:,:,:)                ! needed for zdfmxl 
     821         CALL zdf_mxl( kt, Kmm )                 ! mixed layer depth 
     822         CALL ldf_slp( kt, rhd, rn2, Kbb, Kmm )  ! slopes 
    790823         puslp (:,:,:) = uslp (:,:,:) 
    791824         pvslp (:,:,:) = vslp (:,:,:) 
Note: See TracChangeset for help on using the changeset viewer.