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 12377 for NEMO/trunk/src/OFF/dtadyn.F90 – NEMO

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

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

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

svn merge --ignore-ancestry \

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

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

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

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

Location:
NEMO/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • NEMO/trunk

    • Property svn:externals

      •  

        old new  
        33^/utils/build/mk@HEAD         mk 
        44^/utils/tools@HEAD            tools 
        5 ^/vendors/AGRIF/dev@HEAD      ext/AGRIF 
         5^/vendors/AGRIF/dev_r11615_ENHANCE-04_namelists_as_internalfiles_agrif@HEAD      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 

  • NEMO/trunk/src/OFF/dtadyn.F90

    r11536 r12377  
    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 
    … …  
    4647   PRIVATE 
    4748 
    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 
     49   PUBLIC   dta_dyn_init       ! called by nemo_init 
     50   PUBLIC   dta_dyn            ! called by nemo_gcm 
     51   PUBLIC   dta_dyn_sed_init   ! called by nemo_init 
     52   PUBLIC   dta_dyn_sed        ! called by nemo_gcm 
     53   PUBLIC   dta_dyn_atf        ! called by nemo_gcm 
     54   PUBLIC   dta_dyn_sf_interp  ! called by nemo_gcm 
    5355 
    5456   CHARACTER(len=100) ::   cn_dir          !: Root directory for location of ssr files 
    … …  
    8789   INTEGER, SAVE  :: nprevrec, nsecdyn 
    8890 
     91   !! * Substitutions 
     92#  include "do_loop_substitute.h90" 
    8993   !!---------------------------------------------------------------------- 
    9094   !! NEMO/OFF 4.0 , NEMO Consortium (2018) 
    … …  
    9498CONTAINS 
    9599 
    96    SUBROUTINE dta_dyn( kt ) 
     100   SUBROUTINE dta_dyn( kt, Kbb, Kmm, Kaa ) 
    97101      !!---------------------------------------------------------------------- 
    98102      !!                  ***  ROUTINE dta_dyn  *** 
    … …  
    105109      !!             - interpolates data if needed 
    106110      !!---------------------------------------------------------------------- 
    107       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     111      INTEGER, INTENT(in) ::   kt             ! ocean time-step index 
     112      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa  ! ocean time level indices 
    108113      ! 
    109114      INTEGER             ::   ji, jj, jk 
    … …  
    121126      CALL fld_read( kt, 1, sf_dyn )      !=  read data at kt time step   ==! 
    122127      ! 
    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 
     128      IF( l_ldfslp .AND. .NOT.lk_c1d )   CALL  dta_dyn_slp( kt, Kbb, Kmm )    ! Computation of slopes 
     129      ! 
     130      ts(:,:,:,jp_tem,Kmm) = sf_dyn(jf_tem)%fnow(:,:,:)  * tmask(:,:,:)    ! temperature 
     131      ts(:,:,:,jp_sal,Kmm) = sf_dyn(jf_sal)%fnow(:,:,:)  * tmask(:,:,:)    ! salinity 
    127132      wndm(:,:)         = sf_dyn(jf_wnd)%fnow(:,:,1)  * tmask(:,:,1)    ! wind speed - needed for gas exchange 
    128133      fmmflx(:,:)       = sf_dyn(jf_fmf)%fnow(:,:,1)  * tmask(:,:,1)    ! downward salt flux (v3.5+) 
    … …  
    132137      IF( ln_dynrnf ) THEN  
    133138         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 
     139         IF( ln_dynrnf_depth .AND. .NOT. ln_linssh )    CALL  dta_dyn_hrnf(Kmm) 
     140      ENDIF 
     141      ! 
     142      uu(:,:,:,Kmm)        = sf_dyn(jf_uwd)%fnow(:,:,:) * umask(:,:,:)    ! effective u-transport 
     143      vv(:,:,:,Kmm)        = sf_dyn(jf_vwd)%fnow(:,:,:) * vmask(:,:,:)    ! effective v-transport 
     144      ww(:,:,:)        = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:)    ! effective v-transport 
    140145      ! 
    141146      IF( .NOT.ln_linssh ) THEN 
    … …  
    144149         emp_b  (:,:)   = sf_dyn(jf_empb)%fnow(:,:,1) * tmask(:,:,1)    ! E-P 
    145150         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 
     151         CALL dta_dyn_ssh( kt, zhdivtr, ssh(:,:,Kbb), zemp, ssh(:,:,Kaa), e3t(:,:,:,Kaa) )  !=  ssh, vertical scale factor & vertical transport 
    147152         DEALLOCATE( zemp , zhdivtr ) 
    148153         !                                           Write in the tracer restart file 
    … …  
    152157            IF(lwp) WRITE(numout,*) 'dta_dyn_ssh : ssh field written in tracer restart file at it= ', kt,' date= ', ndastp 
    153158            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    154             CALL iom_rstput( kt, nitrst, numrtw, 'sshn', ssha ) 
    155             CALL iom_rstput( kt, nitrst, numrtw, 'sshb', sshn ) 
     159            CALL iom_rstput( kt, nitrst, numrtw, 'sshn', ssh(:,:,Kaa) ) 
     160            CALL iom_rstput( kt, nitrst, numrtw, 'sshb', ssh(:,:,Kmm) ) 
    156161         ENDIF 
    157162      ENDIF 
    158163      ! 
    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 
     164      CALL eos    ( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 
     165      CALL eos_rab( ts(:,:,:,:,Kmm), rab_n, Kmm )       ! now    local thermal/haline expension ratio at T-points 
     166      CALL bn2    ( ts(:,:,:,:,Kmm), rab_n, rn2, Kmm ) ! before Brunt-Vaisala frequency need for zdfmxl 
     167 
     168      rn2b(:,:,:) = rn2(:,:,:)         ! needed for zdfmxl 
     169      CALL zdf_mxl( kt, Kmm )          ! In any case, we need mxl 
    165170      ! 
    166171      hmld(:,:)       = sf_dyn(jf_mld)%fnow(:,:,1) * tmask(:,:,1)    ! mixed layer depht 
    … …  
    174179      ! 
    175180      ! 
    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   ) 
     181      CALL eos( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 
     182      ! 
     183      IF(sn_cfctl%l_prtctl) THEN                 ! print control 
     184         CALL prt_ctl(tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' tn      - : ', mask1=tmask,  kdim=jpk   ) 
     185         CALL prt_ctl(tab3d_1=ts(:,:,:,jp_sal,Kmm), clinfo1=' sn      - : ', mask1=tmask,  kdim=jpk   ) 
     186         CALL prt_ctl(tab3d_1=uu(:,:,:,Kmm)               , clinfo1=' uu(:,:,:,Kmm)      - : ', mask1=umask,  kdim=jpk   ) 
     187         CALL prt_ctl(tab3d_1=vv(:,:,:,Kmm)               , clinfo1=' vv(:,:,:,Kmm)      - : ', mask1=vmask,  kdim=jpk   ) 
     188         CALL prt_ctl(tab3d_1=ww               , clinfo1=' ww      - : ', mask1=tmask,  kdim=jpk   ) 
    184189         CALL prt_ctl(tab3d_1=avt              , clinfo1=' kz      - : ', mask1=tmask,  kdim=jpk   ) 
    185190         CALL prt_ctl(tab3d_1=uslp             , clinfo1=' slp  - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk) 
    … …  
    192197 
    193198 
    194    SUBROUTINE dta_dyn_init 
     199   SUBROUTINE dta_dyn_init( Kbb, Kmm, Kaa ) 
    195200      !!---------------------------------------------------------------------- 
    196201      !!                  ***  ROUTINE dta_dyn_init  *** 
    … …  
    199204      !! ** Method  : - read the data namdta_dyn namelist 
    200205      !!---------------------------------------------------------------------- 
     206      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa  ! ocean time level indices 
     207      ! 
    201208      INTEGER  :: ierr, ierr0, ierr1, ierr2, ierr3   ! return error code 
    202209      INTEGER  :: ifpr                               ! dummy loop indice 
    … …  
    225232      !!---------------------------------------------------------------------- 
    226233      ! 
    227       REWIND( numnam_ref )              ! Namelist namdta_dyn in reference namelist : Offline: init. of dynamical data 
    228234      READ  ( numnam_ref, namdta_dyn, IOSTAT = ios, ERR = 901) 
    229235901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdta_dyn in reference namelist' ) 
    230       REWIND( numnam_cfg )              ! Namelist namdta_dyn in configuration namelist : Offline: init. of dynamical data 
    231236      READ  ( numnam_cfg, namdta_dyn, IOSTAT = ios, ERR = 902 ) 
    232237902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdta_dyn in configuration namelist' ) 
    … …  
    281286      ! Open file for each variable to get his number of dimension 
    282287      DO ifpr = 1, jfld 
    283          CALL fld_clopn( sf_dyn(ifpr), nyear, nmonth, nday ) 
     288         CALL fld_def( sf_dyn(ifpr) ) 
     289         CALL iom_open( sf_dyn(ifpr)%clname, sf_dyn(ifpr)%num ) 
    284290         idv   = iom_varid( sf_dyn(ifpr)%num , slf_d(ifpr)%clvar )        ! id of the variable sdjf%clvar 
    285291         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 
     292         CALL iom_close( sf_dyn(ifpr)%num )                              ! close file if already open 
    287293         ierr1=0 
    288294         IF( idimv == 3 ) THEN    ! 2D variable 
    … …  
    312318        IF( .NOT. sf_dyn(jf_uwd)%ln_clim .AND. ln_rsttr .AND.    &                     ! Restart: read in restart file 
    313319           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(:,:)   ) 
     320           IF(lwp) WRITE(numout,*) ' ssh(:,:,Kmm) forcing fields read in the restart file for initialisation' 
     321           CALL iom_get( numrtr, jpdom_autoglo, 'sshn', ssh(:,:,Kmm)   ) 
     322           CALL iom_get( numrtr, jpdom_autoglo, 'sshb', ssh(:,:,Kbb)   ) 
    317323        ELSE 
    318            IF(lwp) WRITE(numout,*) ' sshn forcing fields read in the restart file for initialisation' 
     324           IF(lwp) WRITE(numout,*) ' ssh(:,:,Kmm) forcing fields read in the restart file for initialisation' 
    319325           CALL iom_open( 'restart', inum ) 
    320            CALL iom_get( inum, jpdom_autoglo, 'sshn', sshn(:,:)   ) 
    321            CALL iom_get( inum, jpdom_autoglo, 'sshb', sshb(:,:)   ) 
     326           CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh(:,:,Kmm)   ) 
     327           CALL iom_get( inum, jpdom_autoglo, 'sshb', ssh(:,:,Kbb)   ) 
    322328           CALL iom_close( inum )                                        ! close file 
    323329        ENDIF 
    324330        ! 
    325331        DO jk = 1, jpkm1 
    326            e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) ) 
     332           e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( 1._wp + ssh(:,:,Kmm) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) ) 
    327333        ENDDO 
    328         e3t_a(:,:,jpk) = e3t_0(:,:,jpk) 
     334        e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk) 
    329335 
    330336        ! Horizontal scale factor interpolations 
    331337        ! -------------------------------------- 
    332         CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 
    333         CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 
     338        CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 
     339        CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 
    334340 
    335341        ! Vertical scale factor interpolations 
    336342        ! ------------------------------------ 
    337         CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n(:,:,:), 'W' ) 
     343        CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w(:,:,:,Kmm), 'W' ) 
    338344   
    339         e3t_b(:,:,:)  = e3t_n(:,:,:) 
    340         e3u_b(:,:,:)  = e3u_n(:,:,:) 
    341         e3v_b(:,:,:)  = e3v_n(:,:,:) 
     345        e3t(:,:,:,Kbb)  = e3t(:,:,:,Kmm) 
     346        e3u(:,:,:,Kbb)  = e3u(:,:,:,Kmm) 
     347        e3v(:,:,:,Kbb)  = e3v(:,:,:,Kmm) 
    342348 
    343349        ! t- and w- points depth 
    344350        ! ---------------------- 
    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(:,:,:) 
     351        gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 
     352        gdepw(:,:,1,Kmm) = 0.0_wp 
     353 
     354        DO_3D_11_11( 2, jpk ) 
     355          !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere 
     356          !    tmask = wmask, ie everywhere expect at jk = mikt 
     357                                                             ! 1 for jk = 
     358                                                             ! mikt 
     359           zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
     360           gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
     361           gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  & 
     362               &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm)) 
     363        END_3D 
     364 
     365        gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm) 
     366        gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm) 
    365367        ! 
    366368      ENDIF 
    … …  
    374376         ! 
    375377         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 
     378         DO_2D_11_11 
     379            IF( h_rnf(ji,jj) > 0._wp ) THEN 
     380               jk = 2 
     381               DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1 
     382               END DO 
     383               nk_rnf(ji,jj) = jk 
     384            ELSEIF( h_rnf(ji,jj) == -1._wp   ) THEN   ;  nk_rnf(ji,jj) = 1 
     385            ELSEIF( h_rnf(ji,jj) == -999._wp ) THEN   ;  nk_rnf(ji,jj) = mbkt(ji,jj) 
     386            ELSE 
     387               CALL ctl_stop( 'sbc_rnf_init: runoff depth not positive, and not -999 or -1, rnf value in file fort.999'  ) 
     388               WRITE(999,*) 'ji, jj, h_rnf(ji,jj) :', ji, jj, h_rnf(ji,jj) 
     389            ENDIF 
     390         END_2D 
     391         DO_2D_11_11 
     392            h_rnf(ji,jj) = 0._wp 
     393            DO jk = 1, nk_rnf(ji,jj) 
     394               h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm) 
    389395            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 
     396         END_2D 
    399397      ELSE                                       ! runoffs applied at the surface 
    400398         nk_rnf(:,:) = 1 
    401          h_rnf (:,:) = e3t_n(:,:,1) 
     399         h_rnf (:,:) = e3t(:,:,1,Kmm) 
    402400      ENDIF 
    403401      nkrnf_max = MAXVAL( nk_rnf(:,:) ) 
    … …  
    411409      IF(lwp) WRITE(numout,*) ' ' 
    412410      ! 
    413       CALL dta_dyn( nit000 ) 
     411      CALL dta_dyn( nit000, Kbb, Kmm, Kaa ) 
    414412      ! 
    415413   END SUBROUTINE dta_dyn_init 
    416414 
    417    SUBROUTINE dta_dyn_sed( kt ) 
     415   SUBROUTINE dta_dyn_sed( kt, Kmm ) 
    418416      !!---------------------------------------------------------------------- 
    419417      !!                  ***  ROUTINE dta_dyn  *** 
    … …  
    427425      !!---------------------------------------------------------------------- 
    428426      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     427      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index 
    429428      ! 
    430429      !!---------------------------------------------------------------------- 
    … …  
    439438      CALL fld_read( kt, 1, sf_dyn )      !=  read data at kt time step   ==! 
    440439      ! 
    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   ) 
     440      ts(:,:,:,jp_tem,Kmm) = sf_dyn(jf_tem)%fnow(:,:,:)  * tmask(:,:,:)    ! temperature 
     441      ts(:,:,:,jp_sal,Kmm) = sf_dyn(jf_sal)%fnow(:,:,:)  * tmask(:,:,:)    ! salinity 
     442      ! 
     443      CALL eos    ( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 
     444 
     445      IF(sn_cfctl%l_prtctl) THEN                     ! print control 
     446         CALL prt_ctl(tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' tn      - : ', mask1=tmask,  kdim=jpk   ) 
     447         CALL prt_ctl(tab3d_1=ts(:,:,:,jp_sal,Kmm), clinfo1=' sn      - : ', mask1=tmask,  kdim=jpk   ) 
    449448      ENDIF 
    450449      ! 
    … …  
    454453 
    455454 
    456    SUBROUTINE dta_dyn_sed_init 
     455   SUBROUTINE dta_dyn_sed_init( Kmm ) 
    457456      !!---------------------------------------------------------------------- 
    458457      !!                  ***  ROUTINE dta_dyn_init  *** 
    … …  
    461460      !! ** Method  : - read the data namdta_dyn namelist 
    462461      !!---------------------------------------------------------------------- 
     462      INTEGER, INTENT( in ) :: Kmm                   ! ocean time level index 
     463      ! 
    463464      INTEGER  :: ierr, ierr0, ierr1, ierr2, ierr3   ! return error code 
    464465      INTEGER  :: ifpr                               ! dummy loop indice 
    … …  
    475476      !!---------------------------------------------------------------------- 
    476477      ! 
    477       REWIND( numnam_ref )              ! Namelist namdta_dyn in reference namelist : Offline: init. of dynamical data 
    478478      READ  ( numnam_ref, namdta_dyn, IOSTAT = ios, ERR = 901) 
    479479901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdta_dyn in reference namelist' ) 
    480       REWIND( numnam_cfg )              ! Namelist namdta_dyn in configuration namelist : Offline: init. of dynamical data 
    481480      READ  ( numnam_cfg, namdta_dyn, IOSTAT = ios, ERR = 902 ) 
    482481902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdta_dyn in configuration namelist' ) 
    … …  
    508507      ! Open file for each variable to get his number of dimension 
    509508      DO ifpr = 1, jfld 
    510          CALL fld_clopn( sf_dyn(ifpr), nyear, nmonth, nday ) 
     509         CALL fld_def( sf_dyn(ifpr) ) 
     510         CALL iom_open( sf_dyn(ifpr)%clname, sf_dyn(ifpr)%num ) 
    511511         idv   = iom_varid( sf_dyn(ifpr)%num , slf_d(ifpr)%clvar )        ! id of the variable sdjf%clvar 
    512512         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 
     513         CALL iom_close( sf_dyn(ifpr)%num )                              ! close file if already open 
    514514         ierr1=0 
    515515         IF( idimv == 3 ) THEN    ! 2D variable 
    … …  
    525525      END DO 
    526526      ! 
    527       CALL dta_dyn_sed( nit000 ) 
     527      CALL dta_dyn_sed( nit000, Kmm ) 
    528528      ! 
    529529   END SUBROUTINE dta_dyn_sed_init 
    530530 
    531    SUBROUTINE dta_dyn_swp( kt ) 
     531   SUBROUTINE dta_dyn_atf( kt, Kbb, Kmm, Kaa ) 
    532532     !!--------------------------------------------------------------------- 
    533533      !!                    ***  ROUTINE dta_dyn_swp  *** 
    534534      !! 
    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 
     535      !! ** Purpose :   Asselin time filter of now SSH 
     536      !!--------------------------------------------------------------------- 
     537      INTEGER, INTENT(in) :: kt             ! time step 
     538      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa  ! ocean time level indices 
     539      ! 
     540      !!--------------------------------------------------------------------- 
     541 
     542      IF( kt == nit000 ) THEN 
     543         IF(lwp) WRITE(numout,*) 
     544         IF(lwp) WRITE(numout,*) 'dta_dyn_atf : Asselin time filter of sea surface height' 
     545         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 
     546      ENDIF 
     547 
     548      ssh(:,:,Kmm) = ssh(:,:,Kmm) + atfp * ( ssh(:,:,Kbb) - 2 * ssh(:,:,Kmm) + ssh(:,:,Kaa))   
     549 
     550      !! Do we also need to time filter e3t?? 
     551      ! 
     552   END SUBROUTINE dta_dyn_atf 
     553    
     554   SUBROUTINE dta_dyn_sf_interp( kt, Kmm ) 
     555      !!--------------------------------------------------------------------- 
     556      !!                    ***  ROUTINE dta_dyn_sf_interp  *** 
     557      !! 
     558      !! ** Purpose :   Calculate scale factors at U/V/W points and depths 
     559      !!                given the after e3t field 
     560      !!--------------------------------------------------------------------- 
     561      INTEGER, INTENT(in) :: kt   ! time step 
     562      INTEGER, INTENT(in) :: Kmm  ! ocean time level indices 
    539563      ! 
    540564      INTEGER             :: ji, jj, jk 
    … …  
    542566      !!--------------------------------------------------------------------- 
    543567 
    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  
    558568      ! Horizontal scale factor interpolations 
    559569      ! -------------------------------------- 
    560       CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 
    561       CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 
     570      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 
     571      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 
    562572 
    563573      ! Vertical scale factor interpolations 
    564574      ! ------------------------------------ 
    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(:,:,:) 
     575      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' ) 
    570576 
    571577      ! t- and w- points depth 
    572578      ! ---------------------- 
    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     
     579      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 
     580      gdepw(:,:,1,Kmm) = 0.0_wp 
     581      ! 
     582      DO_3D_11_11( 2, jpk ) 
     583         zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
     584         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
     585         gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  & 
     586            &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm)) 
     587      END_3D 
     588      ! 
     589   END SUBROUTINE dta_dyn_sf_interp 
    592590 
    593591   SUBROUTINE dta_dyn_ssh( kt, phdivtr, psshb,  pemp, pssha, pe3ta ) 
    … …  
    595593      !!                ***  ROUTINE dta_dyn_wzv  *** 
    596594      !!                    
    597       !! ** Purpose :   compute the after ssh (ssha) and the now vertical velocity 
     595      !! ** Purpose :   compute the after ssh (ssh(:,:,Kaa)) and the now vertical velocity 
    598596      !! 
    599597      !! ** Method  : Using the incompressibility hypothesis,  
    … …  
    608606      !!          The boundary conditions are w=0 at the bottom (no flux) 
    609607      !! 
    610       !! ** action  :   ssha / e3t_a / wn 
     608      !! ** action  :   ssh(:,:,Kaa) / e3t(:,:,:,Kaa) / ww 
    611609      !! 
    612610      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
    … …  
    641639 
    642640 
    643    SUBROUTINE dta_dyn_hrnf 
     641   SUBROUTINE dta_dyn_hrnf( Kmm ) 
    644642      !!---------------------------------------------------------------------- 
    645643      !!                  ***  ROUTINE sbc_rnf  *** 
    … …  
    654652      !!---------------------------------------------------------------------- 
    655653      !! 
    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 
     654      INTEGER, INTENT(in) :: Kmm  ! ocean time level index 
     655      ! 
     656      INTEGER  ::   ji, jj, jk    ! dummy loop indices 
     657      !!---------------------------------------------------------------------- 
     658      ! 
     659      DO_2D_11_11 
     660         h_rnf(ji,jj) = 0._wp 
     661         DO jk = 1, nk_rnf(ji,jj)                           ! recalculates h_rnf to be the depth in metres 
     662             h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm)   ! to the bottom of the relevant grid box 
     663         END DO 
     664      END_2D 
    667665      ! 
    668666   END SUBROUTINE dta_dyn_hrnf 
    … …  
    670668 
    671669 
    672    SUBROUTINE dta_dyn_slp( kt ) 
     670   SUBROUTINE dta_dyn_slp( kt, Kbb, Kmm ) 
    673671      !!--------------------------------------------------------------------- 
    674672      !!                    ***  ROUTINE dta_dyn_slp  *** 
    … …  
    678676      !!--------------------------------------------------------------------- 
    679677      INTEGER,  INTENT(in) :: kt       ! time step 
     678      INTEGER,  INTENT(in) :: Kbb, Kmm ! ocean time level indices 
    680679      ! 
    681680      INTEGER  ::   ji, jj     ! dummy loop indices 
    … …  
    693692            zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:)   ! salinity  
    694693            avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,1) * tmask(:,:,:)   ! vertical diffusive coef. 
    695             CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
     694            CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm ) 
    696695            uslpdta (:,:,:,1) = zuslp (:,:,:)  
    697696            vslpdta (:,:,:,1) = zvslp (:,:,:)  
    … …  
    702701            zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity  
    703702            avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef. 
    704             CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
     703            CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm ) 
    705704            uslpdta (:,:,:,2) = zuslp (:,:,:)  
    706705            vslpdta (:,:,:,2) = zvslp (:,:,:)  
    … …  
    721720              zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity  
    722721              avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef. 
    723               CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
     722              CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm ) 
    724723              ! 
    725724              uslpdta (:,:,:,2) = zuslp (:,:,:)  
    … …  
    745744         zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:)   ! salinity  
    746745         avt(:,:,:)        = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:)   ! vertical diffusive coef. 
    747          CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
     746         CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm ) 
    748747         ! 
    749748         IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace) 
    … …  
    758757 
    759758 
    760    SUBROUTINE compute_slopes( kt, pts, puslp, pvslp, pwslpi, pwslpj ) 
     759   SUBROUTINE compute_slopes( kt, pts, puslp, pvslp, pwslpi, pwslpj, Kbb, Kmm ) 
    761760      !!--------------------------------------------------------------------- 
    762761      !!                    ***  ROUTINE dta_dyn_slp  *** 
    … …  
    770769      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pwslpi   ! zonal diapycnal slopes 
    771770      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pwslpj   ! meridional diapycnal slopes 
     771      INTEGER ,                              INTENT(in ) :: Kbb, Kmm ! ocean time level indices 
    772772      !!--------------------------------------------------------------------- 
    773773      ! 
    774774      IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace) 
    775775         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 
     776         CALL eos_rab( pts, rab_n, Kmm )       ! now local thermal/haline expension ratio at T-points 
     777         CALL bn2    ( pts, rab_n, rn2, Kmm  ) ! now    Brunt-Vaisala 
    778778 
    779779      ! Partial steps: before Horizontal DErivative 
    780780      IF( ln_zps  .AND. .NOT. ln_isfcav)                            & 
    781          &            CALL zps_hde    ( kt, jpts, pts, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient 
     781         &            CALL zps_hde    ( kt, Kmm, jpts, pts, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient 
    782782         &                                        rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
    783783      IF( ln_zps .AND.        ln_isfcav)                            & 
    784          &            CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv, gtui, gtvi, &  ! Partial steps for top cell (ISF) 
     784         &            CALL zps_hde_isf( kt, Kmm, jpts, pts, gtsu, gtsv, gtui, gtvi, &  ! Partial steps for top cell (ISF) 
    785785         &                                        rhd, gru , grv , grui, grvi )  ! of t, s, rd at the first ocean level 
    786786 
    787          rn2b(:,:,:) = rn2(:,:,:)         ! need for zdfmxl 
    788          CALL zdf_mxl( kt )            ! mixed layer depth 
    789          CALL ldf_slp( kt, rhd, rn2 )  ! slopes 
     787         rn2b(:,:,:) = rn2(:,:,:)                ! needed for zdfmxl 
     788         CALL zdf_mxl( kt, Kmm )                 ! mixed layer depth 
     789         CALL ldf_slp( kt, rhd, rn2, Kbb, Kmm )  ! slopes 
    790790         puslp (:,:,:) = uslp (:,:,:) 
    791791         pvslp (:,:,:) = vslp (:,:,:) 
Note: See TracChangeset for help on using the changeset viewer.