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 – 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:
3 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 (:,:,:) 
  • NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OFF/nemogcm.F90

    r10601 r13463  
    77   !!            3.4  ! 2011-01  (C. Ethe, A. R. Porter, STFC Daresbury) dynamical allocation 
    88   !!            4.0  ! 2016-10  (C. Ethe, G. Madec, S. Flavoni)  domain configuration / user defined interface 
     9   !!            4.1  ! 2019-08  (A. Coward, D. Storkey) rewrite in preparation for new timestepping scheme 
    910   !!---------------------------------------------------------------------- 
    1011 
     
    2728   USE usrdef_nam     ! user defined configuration 
    2829   USE eosbn2         ! equation of state            (eos bn2 routine) 
     30#if defined key_qco 
     31   USE domqco         ! tools for scale factor         (dom_qco_r3c  routine) 
     32#endif 
     33   USE bdyini         ! open boundary cond. setting        (bdy_init routine) 
    2934   !              ! ocean physics 
    3035   USE ldftra         ! lateral diffusivity setting    (ldf_tra_init routine) 
     
    5863   USE timing         ! Timing 
    5964   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
    60    USE lbcnfd  , ONLY : isendto, nsndto, nfsloop, nfeloop   ! Setup of north fold exchanges 
     65   USE lbcnfd  , ONLY : isendto, nsndto   ! Setup of north fold exchanges 
     66   USE step, ONLY : Nbb, Nnn, Naa, Nrhs   ! time level indices 
     67   USE halo_mng 
    6168 
    6269   IMPLICIT NONE 
     
    8895      !!              Madec, 2008, internal report, IPSL. 
    8996      !!---------------------------------------------------------------------- 
    90       INTEGER :: istp, indic       ! time step index 
     97      INTEGER :: istp       ! time step index 
    9198      !!---------------------------------------------------------------------- 
    9299 
     
    111118                                CALL iom_setkt  ( istp - nit000 + 1, cxios_context )   ! say to iom that we are at time step kstp 
    112119#if defined key_sed_off 
    113                                 CALL dta_dyn_sed( istp )         ! Interpolation of the dynamical fields 
     120                                CALL dta_dyn_sed( istp,      Nnn      )       ! Interpolation of the dynamical fields 
    114121#else 
    115                                 CALL dta_dyn    ( istp )         ! Interpolation of the dynamical fields 
    116          IF( .NOT.ln_linssh )   CALL dta_dyn_swp( istp )         ! swap of sea  surface height and vertical scale factors 
    117 #endif 
    118                                 CALL trc_stp    ( istp )         ! time-stepping 
    119                                 CALL stp_ctl    ( istp, indic )  ! Time loop: control and print 
     122                                CALL dta_dyn    ( istp, Nbb, Nnn, Naa )       ! Interpolation of the dynamical fields 
     123#endif 
     124#if ! defined key_sed_off 
     125         IF( .NOT.ln_linssh ) THEN 
     126                                CALL dta_dyn_atf( istp, Nbb, Nnn, Naa )       ! time filter of sea  surface height and vertical scale factors 
     127# if defined key_qco 
     128                                CALL dom_qco_r3c( ssh(:,:,Kmm), r3t_f, r3u_f, r3v_f ) 
     129# endif 
     130         ENDIF 
     131                                CALL trc_stp    ( istp, Nbb, Nnn, Nrhs, Naa ) ! time-stepping 
     132# if defined key_qco 
     133                                !r3t(:,:,Kmm) = r3t_f(:,:)                     ! update ssh to h0 ratio 
     134                                !r3u(:,:,Kmm) = r3u_f(:,:) 
     135                                !r3v(:,:,Kmm) = r3v_f(:,:) 
     136# endif 
     137#endif 
     138         ! Swap time levels 
     139         Nrhs = Nbb 
     140         Nbb = Nnn 
     141         Nnn = Naa 
     142         Naa = Nrhs 
     143         ! 
     144#if ! defined key_qco 
     145#if ! defined key_sed_off 
     146         IF( .NOT.ln_linssh )   CALL dta_dyn_sf_interp( istp, Nnn )  ! calculate now grid parameters 
     147#endif 
     148#endif          
     149                                CALL stp_ctl    ( istp )             ! Time loop: control and print 
    120150         istp = istp + 1 
    121151      END DO 
     
    131161 
    132162      IF( nstop /= 0 .AND. lwp ) THEN                 ! error print 
    133          WRITE(numout,cform_err) 
    134          WRITE(numout,*) '   ==>>>   nemo_gcm: a total of ', nstop, ' errors have been found' 
    135          WRITE(numout,*) 
     163         WRITE(ctmp1,*) '   ==>>>   nemo_gcm: a total of ', nstop, ' errors have been found' 
     164         WRITE(ctmp2,*) '           Look for "E R R O R" messages in all existing ocean_output* files' 
     165         CALL ctl_stop( ' ', ctmp1, ' ', ctmp2 ) 
    136166      ENDIF 
    137167      ! 
     
    146176#endif 
    147177      ! 
     178      IF(lwm) THEN 
     179         IF( nstop == 0 ) THEN   ;   STOP 0 
     180         ELSE                    ;   STOP 123 
     181         ENDIF 
     182      ENDIF 
     183      ! 
    148184   END SUBROUTINE nemo_gcm 
    149185 
     
    155191      !! ** Purpose :   initialization of the nemo model in off-line mode 
    156192      !!---------------------------------------------------------------------- 
    157       INTEGER  ::   ji                 ! dummy loop indices 
    158       INTEGER  ::   ios, ilocal_comm   ! local integers 
    159       CHARACTER(len=120), DIMENSION(30) ::   cltxt, cltxt2, clnam 
    160       !! 
    161       NAMELIST/namctl/ ln_ctl   , sn_cfctl, nn_print, nn_ictls, nn_ictle,   & 
    162          &             nn_isplt , nn_jsplt, nn_jctls, nn_jctle,             & 
    163          &             ln_timing, ln_diacfl 
     193      INTEGER ::   ios, ilocal_comm   ! local integers 
     194      !! 
     195      NAMELIST/namctl/ sn_cfctl, ln_timing, ln_diacfl,                                & 
     196         &             nn_isplt,  nn_jsplt,  nn_ictls, nn_ictle, nn_jctls, nn_jctle 
    164197      NAMELIST/namcfg/ ln_read_cfg, cn_domcfg, ln_closea, ln_write_cfg, cn_domcfg_out, ln_use_jattr 
    165198      !!---------------------------------------------------------------------- 
    166199      ! 
    167       cltxt  = '' 
    168       cltxt2 = '' 
    169       clnam  = ''   
    170200      cxios_context = 'nemo' 
    171       ! 
    172       !                             ! Open reference namelist and configuration namelist files 
    173       CALL ctl_opn( numnam_ref, 'namelist_ref', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE. ) 
    174       CALL ctl_opn( numnam_cfg, 'namelist_cfg', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE. ) 
    175       ! 
    176       REWIND( numnam_ref )              ! Namelist namctl in reference namelist 
     201      nn_hls = 1 
     202      ! 
     203      !                             !-------------------------------------------------! 
     204      !                             !     set communicator & select the local rank    ! 
     205      !                             !  must be done as soon as possible to get narea  ! 
     206      !                             !-------------------------------------------------! 
     207      ! 
     208#if defined key_iomput 
     209      CALL xios_initialize( "for_xios_mpi_id", return_comm=ilocal_comm )   ! nemo local communicator given by xios 
     210      CALL mpp_start( ilocal_comm ) 
     211#else 
     212      CALL mpp_start( ) 
     213#endif 
     214      ! 
     215      narea = mpprank + 1               ! mpprank: the rank of proc (0 --> mppsize -1 ) 
     216      lwm = (narea == 1)                ! control of output namelists 
     217      ! 
     218      !                             !---------------------------------------------------------------! 
     219      !                             ! Open output files, reference and configuration namelist files ! 
     220      !                             !---------------------------------------------------------------! 
     221      ! 
     222      ! open ocean.output as soon as possible to get all output prints (including errors messages) 
     223      IF( lwm )   CALL ctl_opn(     numout,        'ocean.output', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, -1, .FALSE. ) 
     224      ! open reference and configuration namelist files 
     225                  CALL load_nml( numnam_ref,        'namelist_ref',                                           -1, lwm ) 
     226                  CALL load_nml( numnam_cfg,        'namelist_cfg',                                           -1, lwm ) 
     227      IF( lwm )   CALL ctl_opn(     numond, 'output.namelist.dyn', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, -1, .FALSE. ) 
     228      ! open /dev/null file to be able to supress output write easily 
     229      IF( Agrif_Root() ) THEN 
     230                  CALL ctl_opn(     numnul,           '/dev/null', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, -1, .FALSE. ) 
     231#ifdef key_agrif 
     232      ELSE 
     233                  numnul = Agrif_Parent(numnul)    
     234#endif 
     235      ENDIF 
     236      ! 
     237      !                             !--------------------! 
     238      !                             ! Open listing units !  -> need sn_cfctl from namctl to define lwp 
     239      !                             !--------------------! 
     240      ! 
    177241      READ  ( numnam_ref, namctl, IOSTAT = ios, ERR = 901 ) 
    178 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namctl in reference namelist', .TRUE. ) 
    179       REWIND( numnam_cfg )              ! Namelist namctl in confguration namelist 
     242901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namctl in reference namelist' ) 
    180243      READ  ( numnam_cfg, namctl, IOSTAT = ios, ERR = 902 ) 
    181 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namctl in configuration namelist', .TRUE. ) 
    182       ! 
    183       REWIND( numnam_ref )              ! Namelist namcfg in reference namelist 
    184       READ  ( numnam_ref, namcfg, IOSTAT = ios, ERR = 903 ) 
    185 903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namcfg in reference namelist', .TRUE. ) 
    186       REWIND( numnam_cfg )              ! Namelist namcfg in confguration namelist 
    187       READ  ( numnam_cfg, namcfg, IOSTAT = ios, ERR = 904 ) 
    188 904   IF( ios >  0 )   CALL ctl_nam ( ios , 'namcfg in configuration namelist', .TRUE. )    
    189  
    190       !                             !--------------------------! 
    191       !                             !  Set global domain size  !   (control print return in cltxt2) 
    192       !                             !--------------------------! 
    193       IF( ln_read_cfg ) THEN              ! Read sizes in domain configuration file 
    194          CALL domain_cfg ( cltxt2,        cn_cfg, nn_cfg, jpiglo, jpjglo, jpkglo, jperio ) 
    195          ! 
    196       ELSE                                ! user-defined namelist 
    197          CALL usr_def_nam( cltxt2, clnam, cn_cfg, nn_cfg, jpiglo, jpjglo, jpkglo, jperio ) 
    198       ENDIF 
    199       ! 
    200       l_offline = .true.                  ! passive tracers are run offline 
    201       ! 
    202       !                             !--------------------------------------------! 
    203       !                             !  set communicator & select the local node  ! 
    204       !                             !  NB: mynode also opens output.namelist.dyn ! 
    205       !                             !      on unit number numond on first proc   ! 
    206       !                             !--------------------------------------------! 
    207 #if defined key_iomput 
    208       CALL  xios_initialize( "for_xios_mpi_id",return_comm=ilocal_comm ) 
    209       narea = mynode( cltxt, 'output.namelist.dyn', numnam_ref, numnam_cfg, numond , nstop, ilocal_comm )   ! Nodes selection 
    210 #else 
    211       ilocal_comm = 0 
    212       narea = mynode( cltxt, 'output.namelist.dyn', numnam_ref, numnam_cfg, numond , nstop )                ! Nodes selection (control print return in cltxt) 
    213 #endif 
    214  
    215       narea = narea + 1                       ! mynode return the rank of proc (0 --> jpnij -1 ) 
    216  
    217       IF( sn_cfctl%l_config ) THEN 
    218          ! Activate finer control of report outputs 
    219          ! optionally switch off output from selected areas (note this only 
    220          ! applies to output which does not involve global communications) 
    221          IF( ( narea < sn_cfctl%procmin .OR. narea > sn_cfctl%procmax  ) .OR. & 
    222            & ( MOD( narea - sn_cfctl%procmin, sn_cfctl%procincr ) /= 0 ) )    & 
    223            &   CALL nemo_set_cfctl( sn_cfctl, .FALSE., .FALSE. ) 
    224       ELSE 
    225          ! Use ln_ctl to turn on or off all options. 
    226          CALL nemo_set_cfctl( sn_cfctl, ln_ctl, .TRUE. ) 
    227       ENDIF 
    228  
    229       lwm = (narea == 1)                      ! control of output namelists 
    230       lwp = (narea == 1) .OR. ln_ctl          ! control of all listing output print 
    231  
    232       IF(lwm) THEN               ! write merged namelists from earlier to output namelist  
    233          !                       ! now that the file has been opened in call to mynode.  
    234          !                       ! NB: nammpp has already been written in mynode (if lk_mpp_mpi) 
    235          WRITE( numond, namctl ) 
    236          WRITE( numond, namcfg ) 
    237          IF( .NOT.ln_read_cfg ) THEN 
    238             DO ji = 1, SIZE(clnam) 
    239                IF( TRIM(clnam(ji)) /= '' )   WRITE(numond, * ) clnam(ji)     ! namusr_def print 
    240             END DO 
    241          ENDIF 
    242       ENDIF 
    243  
     244902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namctl in configuration namelist' ) 
     245      ! 
     246      ! finalize the definition of namctl variables 
     247      IF( narea < sn_cfctl%procmin .OR. narea > sn_cfctl%procmax .OR. MOD( narea - sn_cfctl%procmin, sn_cfctl%procincr ) /= 0 )   & 
     248         &   CALL nemo_set_cfctl( sn_cfctl, .FALSE. ) 
     249      ! 
     250      lwp = (narea == 1) .OR. sn_cfctl%l_oceout    ! control of all listing output print 
     251      ! 
    244252      IF(lwp) THEN                            ! open listing units 
    245253         ! 
    246          CALL ctl_opn( numout, 'ocean.output', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea ) 
     254         IF( .NOT. lwm )   &           ! alreay opened for narea == 1 
     255            &     CALL ctl_opn(     numout,        'ocean.output', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, -1, .FALSE., narea ) 
    247256         ! 
    248257         WRITE(numout,*) 
    249          WRITE(numout,*) '   CNRS - NERC - Met OFFICE - MERCATOR-ocean - INGV - CMCC' 
     258         WRITE(numout,*) '   CNRS - NERC - Met OFFICE - MERCATOR-ocean - CMCC' 
    250259         WRITE(numout,*) '                       NEMO team' 
    251260         WRITE(numout,*) '                   Off-line TOP Model' 
     
    266275         WRITE(numout,*) "     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ " 
    267276         WRITE(numout,*) 
    268          DO ji = 1, SIZE(cltxt) 
    269             IF( TRIM(cltxt (ji)) /= '' )   WRITE(numout,*) TRIM(cltxt(ji))    ! control print of mynode 
    270          END DO 
    271          WRITE(numout,*) 
    272          WRITE(numout,*) 
    273          DO ji = 1, SIZE(cltxt2) 
    274             IF( TRIM(cltxt2(ji)) /= '' )   WRITE(numout,*) TRIM(cltxt2(ji))   ! control print of domain size 
    275          END DO 
    276277         ! 
    277278         WRITE(numout,cform_aaa)                                        ! Flag AAAAAAA 
    278279         ! 
    279280      ENDIF 
    280       ! open /dev/null file to be able to supress output write easily 
    281       CALL ctl_opn( numnul, '/dev/null', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE. ) 
    282       ! 
    283       !                                      ! Domain decomposition 
    284       CALL mpp_init                          ! MPP 
    285  
     281      ! 
     282      IF(lwm) WRITE( numond, namctl ) 
     283      ! 
     284      !                             !------------------------------------! 
     285      !                             !  Set global domain size parameters ! 
     286      !                             !------------------------------------! 
     287      !      
     288      READ  ( numnam_ref, namcfg, IOSTAT = ios, ERR = 903 ) 
     289903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namcfg in reference namelist' ) 
     290      READ  ( numnam_cfg, namcfg, IOSTAT = ios, ERR = 904 ) 
     291904   IF( ios >  0 )   CALL ctl_nam ( ios , 'namcfg in configuration namelist' )    
     292      ! 
     293      IF( ln_read_cfg ) THEN              ! Read sizes in domain configuration file 
     294         CALL domain_cfg ( cn_cfg, nn_cfg, Ni0glo, Nj0glo, jpkglo, jperio ) 
     295      ELSE                                ! user-defined namelist 
     296         CALL usr_def_nam( cn_cfg, nn_cfg, Ni0glo, Nj0glo, jpkglo, jperio ) 
     297      ENDIF 
     298      ! 
     299      IF(lwm)   WRITE( numond, namcfg ) 
     300      l_offline = .true.                  ! passive tracers are run offline 
     301      ! 
     302      !                             !-----------------------------------------! 
     303      !                             ! mpp parameters and domain decomposition ! 
     304      !                             !-----------------------------------------! 
     305      ! 
     306      CALL mpp_init 
     307 
     308      CALL halo_mng_init() 
    286309      ! Now we know the dimensions of the grid and numout has been set: we can allocate arrays 
    287310      CALL nemo_alloc() 
     311 
     312      ! Initialise time level indices 
     313      Nbb = 1; Nnn = 2; Naa = 3; Nrhs = Naa 
    288314 
    289315      !                             !-------------------------------! 
     
    300326                           CALL     eos_init        ! Equation of state 
    301327      IF( lk_c1d       )   CALL     c1d_init        ! 1D column configuration 
    302                            CALL     dom_init("OPA") ! Domain 
    303       IF( ln_ctl       )   CALL prt_ctl_init        ! Print control 
    304  
    305                            CALL  istate_init    ! ocean initial state (Dynamics and tracers) 
    306  
    307                            CALL     sbc_init    ! Forcings : surface module 
     328                           CALL     dom_init( Nbb, Nnn, Naa, "OPA") ! Domain 
     329      IF( sn_cfctl%l_prtctl )   & 
     330         &                 CALL prt_ctl_init        ! Print control 
     331 
     332                           CALL  istate_init( Nnn, Naa )    ! ocean initial state (Dynamics and tracers) 
     333 
     334                           CALL     sbc_init( Nbb, Nnn, Naa )    ! Forcings : surface module 
     335                           CALL     bdy_init    ! Open boundaries initialisation     
    308336 
    309337      !                                      ! Tracer physics 
     
    319347                           CALL trc_rst_cal( nit000, 'READ' )   ! calendar 
    320348#if defined key_sed_off 
    321                            CALL dta_dyn_sed_init ! Initialization for the dynamics 
     349                           CALL dta_dyn_sed_init(  Nnn      )        ! Initialization for the dynamics 
    322350#else 
    323                            CALL dta_dyn_init   ! Initialization for the dynamics 
    324 #endif 
    325  
    326                            CALL     trc_init   ! Passive tracers initialization 
     351                           CALL dta_dyn_init( Nbb, Nnn, Naa )        ! Initialization for the dynamics 
     352#endif 
     353 
     354                           CALL     trc_init( Nbb, Nnn, Naa )        ! Passive tracers initialization 
    327355                           CALL dia_ptr_init   ! Poleward TRansports initialization 
    328356                            
     
    340368      !! ** Purpose :   control print setting 
    341369      !! 
    342       !! ** Method  : - print namctl information and check some consistencies 
     370      !! ** Method  : - print namctl and namcfg information and check some consistencies 
    343371      !!---------------------------------------------------------------------- 
    344372      ! 
     
    348376         WRITE(numout,*) '~~~~~~~~' 
    349377         WRITE(numout,*) '   Namelist namctl' 
    350          WRITE(numout,*) '      run control (for debugging)     ln_ctl     = ', ln_ctl 
    351          WRITE(numout,*) '       finer control over o/p sn_cfctl%l_config  = ', sn_cfctl%l_config 
    352378         WRITE(numout,*) '                              sn_cfctl%l_runstat = ', sn_cfctl%l_runstat 
    353379         WRITE(numout,*) '                              sn_cfctl%l_trcstat = ', sn_cfctl%l_trcstat 
    354380         WRITE(numout,*) '                              sn_cfctl%l_oceout  = ', sn_cfctl%l_oceout 
    355381         WRITE(numout,*) '                              sn_cfctl%l_layout  = ', sn_cfctl%l_layout 
    356          WRITE(numout,*) '                              sn_cfctl%l_mppout  = ', sn_cfctl%l_mppout 
    357          WRITE(numout,*) '                              sn_cfctl%l_mpptop  = ', sn_cfctl%l_mpptop 
     382         WRITE(numout,*) '                              sn_cfctl%l_prtctl  = ', sn_cfctl%l_prtctl 
     383         WRITE(numout,*) '                              sn_cfctl%l_prttrc  = ', sn_cfctl%l_prttrc 
     384         WRITE(numout,*) '                              sn_cfctl%l_oasout  = ', sn_cfctl%l_oasout 
    358385         WRITE(numout,*) '                              sn_cfctl%procmin   = ', sn_cfctl%procmin   
    359386         WRITE(numout,*) '                              sn_cfctl%procmax   = ', sn_cfctl%procmax   
    360387         WRITE(numout,*) '                              sn_cfctl%procincr  = ', sn_cfctl%procincr  
    361388         WRITE(numout,*) '                              sn_cfctl%ptimincr  = ', sn_cfctl%ptimincr  
    362          WRITE(numout,*) '      level of print                  nn_print   = ', nn_print 
    363          WRITE(numout,*) '      Start i indice for SUM control  nn_ictls   = ', nn_ictls 
    364          WRITE(numout,*) '      End i indice for SUM control    nn_ictle   = ', nn_ictle 
    365          WRITE(numout,*) '      Start j indice for SUM control  nn_jctls   = ', nn_jctls 
    366          WRITE(numout,*) '      End j indice for SUM control    nn_jctle   = ', nn_jctle 
    367          WRITE(numout,*) '      number of proc. following i     nn_isplt   = ', nn_isplt 
    368          WRITE(numout,*) '      number of proc. following j     nn_jsplt   = ', nn_jsplt 
    369389         WRITE(numout,*) '      timing by routine               ln_timing  = ', ln_timing 
    370390         WRITE(numout,*) '      CFL diagnostics                 ln_diacfl  = ', ln_diacfl 
    371391      ENDIF 
    372       ! 
    373       nprint    = nn_print          ! convert DOCTOR namelist names into OLD names 
    374       nictls    = nn_ictls 
    375       nictle    = nn_ictle 
    376       njctls    = nn_jctls 
    377       njctle    = nn_jctle 
    378       isplt     = nn_isplt 
    379       jsplt     = nn_jsplt 
    380  
     392 
     393      IF( .NOT.ln_read_cfg )   ln_closea = .false.   ! dealing possible only with a domcfg file 
    381394      IF(lwp) THEN                  ! control print 
    382395         WRITE(numout,*) 
     
    389402         WRITE(numout,*) '      use file attribute if exists as i/p j-start ln_use_jattr     = ', ln_use_jattr 
    390403      ENDIF 
    391       IF( .NOT.ln_read_cfg )   ln_closea = .false.   ! dealing possible only with a domcfg file 
    392       ! 
    393       !                             ! Parameter control 
    394       ! 
    395       IF( ln_ctl ) THEN                 ! sub-domain area indices for the control prints 
    396          IF( lk_mpp .AND. jpnij > 1 ) THEN 
    397             isplt = jpni   ;   jsplt = jpnj   ;   ijsplt = jpni*jpnj   ! the domain is forced to the real split domain 
    398          ELSE 
    399             IF( isplt == 1 .AND. jsplt == 1  ) THEN 
    400                CALL ctl_warn( ' - isplt & jsplt are equal to 1',   & 
    401                   &           ' - the print control will be done over the whole domain' ) 
    402             ENDIF 
    403             ijsplt = isplt * jsplt            ! total number of processors ijsplt 
    404          ENDIF 
    405          IF(lwp) WRITE(numout,*)'          - The total number of processors over which the' 
    406          IF(lwp) WRITE(numout,*)'            print control will be done is ijsplt : ', ijsplt 
    407          ! 
    408          !                              ! indices used for the SUM control 
    409          IF( nictls+nictle+njctls+njctle == 0 )   THEN    ! print control done over the default area 
    410             lsp_area = .FALSE. 
    411          ELSE                                             ! print control done over a specific  area 
    412             lsp_area = .TRUE. 
    413             IF( nictls < 1 .OR. nictls > jpiglo )   THEN 
    414                CALL ctl_warn( '          - nictls must be 1<=nictls>=jpiglo, it is forced to 1' ) 
    415                nictls = 1 
    416             ENDIF 
    417             IF( nictle < 1 .OR. nictle > jpiglo )   THEN 
    418                CALL ctl_warn( '          - nictle must be 1<=nictle>=jpiglo, it is forced to jpiglo' ) 
    419                nictle = jpiglo 
    420             ENDIF 
    421             IF( njctls < 1 .OR. njctls > jpjglo )   THEN 
    422                CALL ctl_warn( '          - njctls must be 1<=njctls>=jpjglo, it is forced to 1' ) 
    423                njctls = 1 
    424             ENDIF 
    425             IF( njctle < 1 .OR. njctle > jpjglo )   THEN 
    426                CALL ctl_warn( '          - njctle must be 1<=njctle>=jpjglo, it is forced to jpjglo' ) 
    427                njctle = jpjglo 
    428             ENDIF 
    429          ENDIF 
    430       ENDIF 
    431404      ! 
    432405      IF( 1._wp /= SIGN(1._wp,-0._wp)  )   CALL ctl_stop( 'nemo_ctl: The intrinsec SIGN function follows f2003 standard.',  & 
     
    448421      ! 
    449422      IF( numstp     /= -1 )   CLOSE( numstp     )   ! time-step file 
    450       IF( numnam_ref /= -1 )   CLOSE( numnam_ref )   ! oce reference namelist 
    451       IF( numnam_cfg /= -1 )   CLOSE( numnam_cfg )   ! oce configuration namelist 
    452       IF( numout     /=  6 )   CLOSE( numout     )   ! standard model output file 
    453423      IF( lwm.AND.numond  /= -1 )   CLOSE( numond          )   ! oce output namelist 
    454424      ! 
     
    470440      USE zdf_oce,   ONLY : zdf_oce_alloc 
    471441      USE trc_oce,   ONLY : trc_oce_alloc 
     442      USE bdy_oce,   ONLY : bdy_oce_alloc 
    472443      ! 
    473444      INTEGER :: ierr 
     
    479450      ierr = ierr + zdf_oce_alloc()          ! ocean vertical physics 
    480451      ierr = ierr + trc_oce_alloc()          ! shared TRC / TRA arrays 
     452      ierr = ierr + bdy_oce_alloc()    ! bdy masks (incl. initialization)       
    481453      ! 
    482454      CALL mpp_sum( 'nemogcm', ierr ) 
     
    485457   END SUBROUTINE nemo_alloc 
    486458 
    487    SUBROUTINE nemo_set_cfctl(sn_cfctl, setto, for_all ) 
     459   SUBROUTINE nemo_set_cfctl(sn_cfctl, setto ) 
    488460      !!---------------------------------------------------------------------- 
    489461      !!                     ***  ROUTINE nemo_set_cfctl  *** 
    490462      !! 
    491463      !! ** Purpose :   Set elements of the output control structure to setto. 
    492       !!                for_all should be .false. unless all areas are to be 
    493       !!                treated identically. 
    494       !! 
     464     !! 
    495465      !! ** Method  :   Note this routine can be used to switch on/off some 
    496       !!                types of output for selected areas but any output types 
    497       !!                that involve global communications (e.g. mpp_max, glob_sum) 
    498       !!                should be protected from selective switching by the 
    499       !!                for_all argument 
    500       !!---------------------------------------------------------------------- 
    501       LOGICAL :: setto, for_all 
    502       TYPE(sn_ctl) :: sn_cfctl 
    503       !!---------------------------------------------------------------------- 
    504       IF( for_all ) THEN 
    505          sn_cfctl%l_runstat = setto 
    506          sn_cfctl%l_trcstat = setto 
    507       ENDIF 
     466      !!                types of output for selected areas. 
     467      !!---------------------------------------------------------------------- 
     468      TYPE(sn_ctl), INTENT(inout) :: sn_cfctl 
     469      LOGICAL     , INTENT(in   ) :: setto 
     470      !!---------------------------------------------------------------------- 
     471      sn_cfctl%l_runstat = setto 
     472      sn_cfctl%l_trcstat = setto 
    508473      sn_cfctl%l_oceout  = setto 
    509474      sn_cfctl%l_layout  = setto 
    510       sn_cfctl%l_mppout  = setto 
    511       sn_cfctl%l_mpptop  = setto 
     475      sn_cfctl%l_prtctl  = setto 
     476      sn_cfctl%l_prttrc  = setto 
     477      sn_cfctl%l_oasout  = setto 
    512478   END SUBROUTINE nemo_set_cfctl 
    513479 
    514    SUBROUTINE istate_init 
     480   SUBROUTINE istate_init( Kmm, Kaa ) 
    515481      !!---------------------------------------------------------------------- 
    516482      !!                   ***  ROUTINE istate_init  *** 
     
    518484      !! ** Purpose :   Initialization to zero of the dynamics and tracers. 
    519485      !!---------------------------------------------------------------------- 
     486      INTEGER, INTENT(in) ::   Kmm, Kaa  ! ocean time level indices 
    520487      ! 
    521488      !     now fields         !     after fields      ! 
    522       un   (:,:,:)   = 0._wp   ;   ua(:,:,:) = 0._wp   ! 
    523       vn   (:,:,:)   = 0._wp   ;   va(:,:,:) = 0._wp   ! 
    524       wn   (:,:,:)   = 0._wp   !                       ! 
    525       hdivn(:,:,:)   = 0._wp   !                       ! 
    526       tsn  (:,:,:,:) = 0._wp   !                       ! 
     489      uu   (:,:,:,Kmm)   = 0._wp   ;   uu(:,:,:,Kaa) = 0._wp   ! 
     490      vv   (:,:,:,Kmm)   = 0._wp   ;   vv(:,:,:,Kaa) = 0._wp   ! 
     491      ww   (:,:,:)   = 0._wp   !                       ! 
     492      hdiv (:,:,:)   = 0._wp   !                       ! 
     493      ts  (:,:,:,:,Kmm) = 0._wp   !                       ! 
    527494      ! 
    528495      rhd  (:,:,:) = 0.e0 
     
    533500 
    534501 
    535    SUBROUTINE stp_ctl( kt, kindic ) 
     502   SUBROUTINE stp_ctl( kt ) 
    536503      !!---------------------------------------------------------------------- 
    537504      !!                    ***  ROUTINE stp_ctl  *** 
     
    544511      !!---------------------------------------------------------------------- 
    545512      INTEGER, INTENT(in   ) ::   kt      ! ocean time-step index 
    546       INTEGER, INTENT(inout) ::   kindic  ! indicator of solver convergence 
    547513      !!---------------------------------------------------------------------- 
    548514      ! 
Note: See TracChangeset for help on using the changeset viewer.