Ignore:
Timestamp:
2019-11-22T15:29:17+01:00 (19 months ago)
Author:
acc
Message:

Merge in changes from 2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps. This just creates a fresh copy of this branch to use as the merge base. See ticket #2341

Location:
NEMO/branches/2019/dev_r11943_MERGE_2019/src
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11943_MERGE_2019/src

    • Property svn:mergeinfo deleted
  • NEMO/branches/2019/dev_r11943_MERGE_2019/src/OFF/dtadyn.F90

    r11536 r11949  
    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 
     
    9496CONTAINS 
    9597 
    96    SUBROUTINE dta_dyn( kt ) 
     98   SUBROUTINE dta_dyn( kt, Kbb, Kmm, Kaa ) 
    9799      !!---------------------------------------------------------------------- 
    98100      !!                  ***  ROUTINE dta_dyn  *** 
     
    105107      !!             - interpolates data if needed 
    106108      !!---------------------------------------------------------------------- 
    107       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     109      INTEGER, INTENT(in) ::   kt             ! ocean time-step index 
     110      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa  ! ocean time level indices 
    108111      ! 
    109112      INTEGER             ::   ji, jj, jk 
     
    121124      CALL fld_read( kt, 1, sf_dyn )      !=  read data at kt time step   ==! 
    122125      ! 
    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 
     126      IF( l_ldfslp .AND. .NOT.lk_c1d )   CALL  dta_dyn_slp( kt, Kbb, Kmm )    ! Computation of slopes 
     127      ! 
     128      ts(:,:,:,jp_tem,Kmm) = sf_dyn(jf_tem)%fnow(:,:,:)  * tmask(:,:,:)    ! temperature 
     129      ts(:,:,:,jp_sal,Kmm) = sf_dyn(jf_sal)%fnow(:,:,:)  * tmask(:,:,:)    ! salinity 
    127130      wndm(:,:)         = sf_dyn(jf_wnd)%fnow(:,:,1)  * tmask(:,:,1)    ! wind speed - needed for gas exchange 
    128131      fmmflx(:,:)       = sf_dyn(jf_fmf)%fnow(:,:,1)  * tmask(:,:,1)    ! downward salt flux (v3.5+) 
     
    132135      IF( ln_dynrnf ) THEN  
    133136         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 
     137         IF( ln_dynrnf_depth .AND. .NOT. ln_linssh )    CALL  dta_dyn_hrnf(Kmm) 
     138      ENDIF 
     139      ! 
     140      uu(:,:,:,Kmm)        = sf_dyn(jf_uwd)%fnow(:,:,:) * umask(:,:,:)    ! effective u-transport 
     141      vv(:,:,:,Kmm)        = sf_dyn(jf_vwd)%fnow(:,:,:) * vmask(:,:,:)    ! effective v-transport 
     142      ww(:,:,:)        = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:)    ! effective v-transport 
    140143      ! 
    141144      IF( .NOT.ln_linssh ) THEN 
     
    144147         emp_b  (:,:)   = sf_dyn(jf_empb)%fnow(:,:,1) * tmask(:,:,1)    ! E-P 
    145148         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 
     149         CALL dta_dyn_ssh( kt, zhdivtr, ssh(:,:,Kbb), zemp, ssh(:,:,Kaa), e3t(:,:,:,Kaa) )  !=  ssh, vertical scale factor & vertical transport 
    147150         DEALLOCATE( zemp , zhdivtr ) 
    148151         !                                           Write in the tracer restart file 
     
    152155            IF(lwp) WRITE(numout,*) 'dta_dyn_ssh : ssh field written in tracer restart file at it= ', kt,' date= ', ndastp 
    153156            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    154             CALL iom_rstput( kt, nitrst, numrtw, 'sshn', ssha ) 
    155             CALL iom_rstput( kt, nitrst, numrtw, 'sshb', sshn ) 
     157            CALL iom_rstput( kt, nitrst, numrtw, 'sshn', ssh(:,:,Kaa) ) 
     158            CALL iom_rstput( kt, nitrst, numrtw, 'sshb', ssh(:,:,Kmm) ) 
    156159         ENDIF 
    157160      ENDIF 
    158161      ! 
    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 
     162      CALL eos    ( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 
     163      CALL eos_rab( ts(:,:,:,:,Kmm), rab_n, Kmm )       ! now    local thermal/haline expension ratio at T-points 
     164      CALL bn2    ( ts(:,:,:,:,Kmm), rab_n, rn2, Kmm ) ! before Brunt-Vaisala frequency need for zdfmxl 
     165 
     166      rn2b(:,:,:) = rn2(:,:,:)         ! needed for zdfmxl 
     167      CALL zdf_mxl( kt, Kmm )          ! In any case, we need mxl 
    165168      ! 
    166169      hmld(:,:)       = sf_dyn(jf_mld)%fnow(:,:,1) * tmask(:,:,1)    ! mixed layer depht 
     
    174177      ! 
    175178      ! 
    176       CALL eos( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 
     179      CALL eos( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 
    177180      ! 
    178181      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   ) 
     182         CALL prt_ctl(tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' tn      - : ', mask1=tmask,  kdim=jpk   ) 
     183         CALL prt_ctl(tab3d_1=ts(:,:,:,jp_sal,Kmm), clinfo1=' sn      - : ', mask1=tmask,  kdim=jpk   ) 
     184         CALL prt_ctl(tab3d_1=uu(:,:,:,Kmm)               , clinfo1=' uu(:,:,:,Kmm)      - : ', mask1=umask,  kdim=jpk   ) 
     185         CALL prt_ctl(tab3d_1=vv(:,:,:,Kmm)               , clinfo1=' vv(:,:,:,Kmm)      - : ', mask1=vmask,  kdim=jpk   ) 
     186         CALL prt_ctl(tab3d_1=ww               , clinfo1=' ww      - : ', mask1=tmask,  kdim=jpk   ) 
    184187         CALL prt_ctl(tab3d_1=avt              , clinfo1=' kz      - : ', mask1=tmask,  kdim=jpk   ) 
    185188         CALL prt_ctl(tab3d_1=uslp             , clinfo1=' slp  - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk) 
     
    192195 
    193196 
    194    SUBROUTINE dta_dyn_init 
     197   SUBROUTINE dta_dyn_init( Kbb, Kmm, Kaa ) 
    195198      !!---------------------------------------------------------------------- 
    196199      !!                  ***  ROUTINE dta_dyn_init  *** 
     
    199202      !! ** Method  : - read the data namdta_dyn namelist 
    200203      !!---------------------------------------------------------------------- 
     204      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa  ! ocean time level indices 
     205      ! 
    201206      INTEGER  :: ierr, ierr0, ierr1, ierr2, ierr3   ! return error code 
    202207      INTEGER  :: ifpr                               ! dummy loop indice 
     
    312317        IF( .NOT. sf_dyn(jf_uwd)%ln_clim .AND. ln_rsttr .AND.    &                     ! Restart: read in restart file 
    313318           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(:,:)   ) 
     319           IF(lwp) WRITE(numout,*) ' ssh(:,:,Kmm) forcing fields read in the restart file for initialisation' 
     320           CALL iom_get( numrtr, jpdom_autoglo, 'sshn', ssh(:,:,Kmm)   ) 
     321           CALL iom_get( numrtr, jpdom_autoglo, 'sshb', ssh(:,:,Kbb)   ) 
    317322        ELSE 
    318            IF(lwp) WRITE(numout,*) ' sshn forcing fields read in the restart file for initialisation' 
     323           IF(lwp) WRITE(numout,*) ' ssh(:,:,Kmm) forcing fields read in the restart file for initialisation' 
    319324           CALL iom_open( 'restart', inum ) 
    320            CALL iom_get( inum, jpdom_autoglo, 'sshn', sshn(:,:)   ) 
    321            CALL iom_get( inum, jpdom_autoglo, 'sshb', sshb(:,:)   ) 
     325           CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh(:,:,Kmm)   ) 
     326           CALL iom_get( inum, jpdom_autoglo, 'sshb', ssh(:,:,Kbb)   ) 
    322327           CALL iom_close( inum )                                        ! close file 
    323328        ENDIF 
    324329        ! 
    325330        DO jk = 1, jpkm1 
    326            e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) ) 
     331           e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( 1._wp + ssh(:,:,Kmm) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) ) 
    327332        ENDDO 
    328         e3t_a(:,:,jpk) = e3t_0(:,:,jpk) 
     333        e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk) 
    329334 
    330335        ! Horizontal scale factor interpolations 
    331336        ! -------------------------------------- 
    332         CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 
    333         CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 
     337        CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 
     338        CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 
    334339 
    335340        ! Vertical scale factor interpolations 
    336341        ! ------------------------------------ 
    337         CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n(:,:,:), 'W' ) 
     342        CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w(:,:,:,Kmm), 'W' ) 
    338343   
    339         e3t_b(:,:,:)  = e3t_n(:,:,:) 
    340         e3u_b(:,:,:)  = e3u_n(:,:,:) 
    341         e3v_b(:,:,:)  = e3v_n(:,:,:) 
     344        e3t(:,:,:,Kbb)  = e3t(:,:,:,Kmm) 
     345        e3u(:,:,:,Kbb)  = e3u(:,:,:,Kmm) 
     346        e3v(:,:,:,Kbb)  = e3v(:,:,:,Kmm) 
    342347 
    343348        ! t- and w- points depth 
    344349        ! ---------------------- 
    345         gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 
    346         gdepw_n(:,:,1) = 0.0_wp 
     350        gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 
     351        gdepw(:,:,1,Kmm) = 0.0_wp 
    347352 
    348353        DO jk = 2, jpk 
     
    354359                                                                   ! mikt 
    355360                 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)) 
     361                 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
     362                 gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  & 
     363                     &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm)) 
    359364              END DO 
    360365           END DO 
    361366        END DO 
    362367 
    363         gdept_b(:,:,:) = gdept_n(:,:,:) 
    364         gdepw_b(:,:,:) = gdepw_n(:,:,:) 
     368        gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm) 
     369        gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm) 
    365370        ! 
    366371      ENDIF 
     
    393398               h_rnf(ji,jj) = 0._wp 
    394399               DO jk = 1, nk_rnf(ji,jj) 
    395                   h_rnf(ji,jj) = h_rnf(ji,jj) + e3t_n(ji,jj,jk) 
     400                  h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm) 
    396401               END DO 
    397402            END DO 
     
    399404      ELSE                                       ! runoffs applied at the surface 
    400405         nk_rnf(:,:) = 1 
    401          h_rnf (:,:) = e3t_n(:,:,1) 
     406         h_rnf (:,:) = e3t(:,:,1,Kmm) 
    402407      ENDIF 
    403408      nkrnf_max = MAXVAL( nk_rnf(:,:) ) 
     
    411416      IF(lwp) WRITE(numout,*) ' ' 
    412417      ! 
    413       CALL dta_dyn( nit000 ) 
     418      CALL dta_dyn( nit000, Kbb, Kmm, Kaa ) 
    414419      ! 
    415420   END SUBROUTINE dta_dyn_init 
    416421 
    417    SUBROUTINE dta_dyn_sed( kt ) 
     422   SUBROUTINE dta_dyn_sed( kt, Kmm ) 
    418423      !!---------------------------------------------------------------------- 
    419424      !!                  ***  ROUTINE dta_dyn  *** 
     
    427432      !!---------------------------------------------------------------------- 
    428433      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     434      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index 
    429435      ! 
    430436      !!---------------------------------------------------------------------- 
     
    439445      CALL fld_read( kt, 1, sf_dyn )      !=  read data at kt time step   ==! 
    440446      ! 
    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 
     447      ts(:,:,:,jp_tem,Kmm) = sf_dyn(jf_tem)%fnow(:,:,:)  * tmask(:,:,:)    ! temperature 
     448      ts(:,:,:,jp_sal,Kmm) = sf_dyn(jf_sal)%fnow(:,:,:)  * tmask(:,:,:)    ! salinity 
     449      ! 
     450      CALL eos    ( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 
    445451 
    446452      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   ) 
     453         CALL prt_ctl(tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' tn      - : ', mask1=tmask,  kdim=jpk   ) 
     454         CALL prt_ctl(tab3d_1=ts(:,:,:,jp_sal,Kmm), clinfo1=' sn      - : ', mask1=tmask,  kdim=jpk   ) 
    449455      ENDIF 
    450456      ! 
     
    454460 
    455461 
    456    SUBROUTINE dta_dyn_sed_init 
     462   SUBROUTINE dta_dyn_sed_init( Kmm ) 
    457463      !!---------------------------------------------------------------------- 
    458464      !!                  ***  ROUTINE dta_dyn_init  *** 
     
    461467      !! ** Method  : - read the data namdta_dyn namelist 
    462468      !!---------------------------------------------------------------------- 
     469      INTEGER, INTENT( in ) :: Kmm                   ! ocean time level index 
     470      ! 
    463471      INTEGER  :: ierr, ierr0, ierr1, ierr2, ierr3   ! return error code 
    464472      INTEGER  :: ifpr                               ! dummy loop indice 
     
    525533      END DO 
    526534      ! 
    527       CALL dta_dyn_sed( nit000 ) 
     535      CALL dta_dyn_sed( nit000, Kmm ) 
    528536      ! 
    529537   END SUBROUTINE dta_dyn_sed_init 
    530538 
    531    SUBROUTINE dta_dyn_swp( kt ) 
     539   SUBROUTINE dta_dyn_atf( kt, Kbb, Kmm, Kaa ) 
    532540     !!--------------------------------------------------------------------- 
    533541      !!                    ***  ROUTINE dta_dyn_swp  *** 
    534542      !! 
    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 
     543      !! ** Purpose :   Asselin time filter of now SSH 
     544      !!--------------------------------------------------------------------- 
     545      INTEGER, INTENT(in) :: kt             ! time step 
     546      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa  ! ocean time level indices 
     547      ! 
     548      !!--------------------------------------------------------------------- 
     549 
     550      IF( kt == nit000 ) THEN 
     551         IF(lwp) WRITE(numout,*) 
     552         IF(lwp) WRITE(numout,*) 'dta_dyn_atf : Asselin time filter of sea surface height' 
     553         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 
     554      ENDIF 
     555 
     556      ssh(:,:,Kmm) = ssh(:,:,Kmm) + atfp * ( ssh(:,:,Kbb) - 2 * ssh(:,:,Kmm) + ssh(:,:,Kaa))   
     557 
     558      !! Do we also need to time filter e3t?? 
     559      ! 
     560   END SUBROUTINE dta_dyn_atf 
     561    
     562   SUBROUTINE dta_dyn_sf_interp( kt, Kmm ) 
     563      !!--------------------------------------------------------------------- 
     564      !!                    ***  ROUTINE dta_dyn_sf_interp  *** 
     565      !! 
     566      !! ** Purpose :   Calculate scale factors at U/V/W points and depths 
     567      !!                given the after e3t field 
     568      !!--------------------------------------------------------------------- 
     569      INTEGER, INTENT(in) :: kt   ! time step 
     570      INTEGER, INTENT(in) :: Kmm  ! ocean time level indices 
    539571      ! 
    540572      INTEGER             :: ji, jj, jk 
     
    542574      !!--------------------------------------------------------------------- 
    543575 
    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  
    558576      ! Horizontal scale factor interpolations 
    559577      ! -------------------------------------- 
    560       CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 
    561       CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 
     578      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 
     579      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 
    562580 
    563581      ! Vertical scale factor interpolations 
    564582      ! ------------------------------------ 
    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(:,:,:) 
     583      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' ) 
    570584 
    571585      ! t- and w- points depth 
    572586      ! ---------------------- 
    573       gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 
    574       gdepw_n(:,:,1) = 0.0_wp 
     587      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 
     588      gdepw(:,:,1,Kmm) = 0.0_wp 
    575589      ! 
    576590      DO jk = 2, jpk 
     
    578592            DO ji = 1,jpi 
    579593               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)) 
     594               gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
     595               gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  & 
     596                  &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm)) 
    583597            END DO 
    584598         END DO 
    585599      END DO 
    586600      ! 
    587       gdept_b(:,:,:) = gdept_n(:,:,:) 
    588       gdepw_b(:,:,:) = gdepw_n(:,:,:) 
    589       ! 
    590    END SUBROUTINE dta_dyn_swp 
    591     
     601   END SUBROUTINE dta_dyn_sf_interp 
    592602 
    593603   SUBROUTINE dta_dyn_ssh( kt, phdivtr, psshb,  pemp, pssha, pe3ta ) 
     
    595605      !!                ***  ROUTINE dta_dyn_wzv  *** 
    596606      !!                    
    597       !! ** Purpose :   compute the after ssh (ssha) and the now vertical velocity 
     607      !! ** Purpose :   compute the after ssh (ssh(:,:,Kaa)) and the now vertical velocity 
    598608      !! 
    599609      !! ** Method  : Using the incompressibility hypothesis,  
     
    608618      !!          The boundary conditions are w=0 at the bottom (no flux) 
    609619      !! 
    610       !! ** action  :   ssha / e3t_a / wn 
     620      !! ** action  :   ssh(:,:,Kaa) / e3t(:,:,:,Kaa) / ww 
    611621      !! 
    612622      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
     
    641651 
    642652 
    643    SUBROUTINE dta_dyn_hrnf 
     653   SUBROUTINE dta_dyn_hrnf( Kmm ) 
    644654      !!---------------------------------------------------------------------- 
    645655      !!                  ***  ROUTINE sbc_rnf  *** 
     
    654664      !!---------------------------------------------------------------------- 
    655665      !! 
    656       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     666      INTEGER, INTENT(in) :: Kmm  ! ocean time level index 
     667      ! 
     668      INTEGER  ::   ji, jj, jk    ! dummy loop indices 
    657669      !!---------------------------------------------------------------------- 
    658670      ! 
     
    661673            h_rnf(ji,jj) = 0._wp 
    662674            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 
     675                h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm)   ! to the bottom of the relevant grid box 
    664676            END DO 
    665677        END DO 
     
    670682 
    671683 
    672    SUBROUTINE dta_dyn_slp( kt ) 
     684   SUBROUTINE dta_dyn_slp( kt, Kbb, Kmm ) 
    673685      !!--------------------------------------------------------------------- 
    674686      !!                    ***  ROUTINE dta_dyn_slp  *** 
     
    678690      !!--------------------------------------------------------------------- 
    679691      INTEGER,  INTENT(in) :: kt       ! time step 
     692      INTEGER,  INTENT(in) :: Kbb, Kmm ! ocean time level indices 
    680693      ! 
    681694      INTEGER  ::   ji, jj     ! dummy loop indices 
     
    693706            zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:)   ! salinity  
    694707            avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,1) * tmask(:,:,:)   ! vertical diffusive coef. 
    695             CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
     708            CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm ) 
    696709            uslpdta (:,:,:,1) = zuslp (:,:,:)  
    697710            vslpdta (:,:,:,1) = zvslp (:,:,:)  
     
    702715            zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity  
    703716            avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef. 
    704             CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
     717            CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm ) 
    705718            uslpdta (:,:,:,2) = zuslp (:,:,:)  
    706719            vslpdta (:,:,:,2) = zvslp (:,:,:)  
     
    721734              zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity  
    722735              avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef. 
    723               CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
     736              CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm ) 
    724737              ! 
    725738              uslpdta (:,:,:,2) = zuslp (:,:,:)  
     
    745758         zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:)   ! salinity  
    746759         avt(:,:,:)        = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:)   ! vertical diffusive coef. 
    747          CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 
     760         CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj, Kbb, Kmm ) 
    748761         ! 
    749762         IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace) 
     
    758771 
    759772 
    760    SUBROUTINE compute_slopes( kt, pts, puslp, pvslp, pwslpi, pwslpj ) 
     773   SUBROUTINE compute_slopes( kt, pts, puslp, pvslp, pwslpi, pwslpj, Kbb, Kmm ) 
    761774      !!--------------------------------------------------------------------- 
    762775      !!                    ***  ROUTINE dta_dyn_slp  *** 
     
    770783      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pwslpi   ! zonal diapycnal slopes 
    771784      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(out) :: pwslpj   ! meridional diapycnal slopes 
     785      INTEGER ,                              INTENT(in ) :: Kbb, Kmm ! ocean time level indices 
    772786      !!--------------------------------------------------------------------- 
    773787      ! 
    774788      IF( l_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace) 
    775789         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 
     790         CALL eos_rab( pts, rab_n, Kmm )       ! now local thermal/haline expension ratio at T-points 
     791         CALL bn2    ( pts, rab_n, rn2, Kmm  ) ! now    Brunt-Vaisala 
    778792 
    779793      ! Partial steps: before Horizontal DErivative 
    780794      IF( ln_zps  .AND. .NOT. ln_isfcav)                            & 
    781          &            CALL zps_hde    ( kt, jpts, pts, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient 
     795         &            CALL zps_hde    ( kt, Kmm, jpts, pts, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient 
    782796         &                                        rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
    783797      IF( ln_zps .AND.        ln_isfcav)                            & 
    784          &            CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv, gtui, gtvi, &  ! Partial steps for top cell (ISF) 
     798         &            CALL zps_hde_isf( kt, Kmm, jpts, pts, gtsu, gtsv, gtui, gtvi, &  ! Partial steps for top cell (ISF) 
    785799         &                                        rhd, gru , grv , grui, grvi )  ! of t, s, rd at the first ocean level 
    786800 
    787          rn2b(:,:,:) = rn2(:,:,:)         ! need for zdfmxl 
    788          CALL zdf_mxl( kt )            ! mixed layer depth 
    789          CALL ldf_slp( kt, rhd, rn2 )  ! slopes 
     801         rn2b(:,:,:) = rn2(:,:,:)                ! needed for zdfmxl 
     802         CALL zdf_mxl( kt, Kmm )                 ! mixed layer depth 
     803         CALL ldf_slp( kt, rhd, rn2, Kbb, Kmm )  ! slopes 
    790804         puslp (:,:,:) = uslp (:,:,:) 
    791805         pvslp (:,:,:) = vslp (:,:,:) 
  • NEMO/branches/2019/dev_r11943_MERGE_2019/src/OFF/nemogcm.F90

    r11536 r11949  
    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 
     
    5960   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
    6061   USE lbcnfd  , ONLY : isendto, nsndto, nfsloop, nfeloop   ! Setup of north fold exchanges 
     62   USE step, ONLY : Nbb, Nnn, Naa, Nrhs   ! time level indices 
    6163 
    6264   IMPLICIT NONE 
     
    111113                                CALL iom_setkt  ( istp - nit000 + 1, cxios_context )   ! say to iom that we are at time step kstp 
    112114#if defined key_sed_off 
    113                                 CALL dta_dyn_sed( istp )         ! Interpolation of the dynamical fields 
     115                                CALL dta_dyn_sed( istp,      Nnn      )       ! Interpolation of the dynamical fields 
    114116#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 
     117                                CALL dta_dyn    ( istp, Nbb, Nnn, Naa )       ! Interpolation of the dynamical fields 
     118#endif 
     119                                CALL trc_stp    ( istp, Nbb, Nnn, Nrhs, Naa ) ! time-stepping 
     120#if ! defined key_sed_off 
     121         IF( .NOT.ln_linssh )   CALL dta_dyn_atf( istp, Nbb, Nnn, Naa )       ! time filter of sea  surface height and vertical scale factors 
     122#endif 
     123         ! Swap time levels 
     124         Nrhs = Nbb 
     125         Nbb = Nnn 
     126         Nnn = Naa 
     127         Naa = Nrhs 
     128         ! 
     129#if ! defined key_sed_off 
     130         IF( .NOT.ln_linssh )   CALL dta_dyn_sf_interp( istp, Nnn )  ! calculate now grid parameters 
     131#endif 
    119132                                CALL stp_ctl    ( istp, indic )  ! Time loop: control and print 
    120133         istp = istp + 1 
     
    285298      CALL nemo_alloc() 
    286299 
     300      ! Initialise time level indices 
     301      Nbb = 1; Nnn = 2; Naa = 3; Nrhs = Naa 
     302    
     303 
    287304      !                             !-------------------------------! 
    288305      !                             !  NEMO general initialization  ! 
     
    298315                           CALL     eos_init        ! Equation of state 
    299316      IF( lk_c1d       )   CALL     c1d_init        ! 1D column configuration 
    300                            CALL     dom_init("OPA") ! Domain 
     317                           CALL     dom_init( Nbb, Nnn, Naa, "OPA") ! Domain 
    301318      IF( ln_ctl       )   CALL prt_ctl_init        ! Print control 
    302319 
    303                            CALL  istate_init    ! ocean initial state (Dynamics and tracers) 
    304  
    305                            CALL     sbc_init    ! Forcings : surface module 
     320                           CALL  istate_init( Nnn, Naa )    ! ocean initial state (Dynamics and tracers) 
     321 
     322                           CALL     sbc_init( Nbb, Nnn, Naa )    ! Forcings : surface module 
    306323 
    307324      !                                      ! Tracer physics 
     
    317334                           CALL trc_rst_cal( nit000, 'READ' )   ! calendar 
    318335#if defined key_sed_off 
    319                            CALL dta_dyn_sed_init ! Initialization for the dynamics 
     336                           CALL dta_dyn_sed_init(  Nnn      )        ! Initialization for the dynamics 
    320337#else 
    321                            CALL dta_dyn_init   ! Initialization for the dynamics 
    322 #endif 
    323  
    324                            CALL     trc_init   ! Passive tracers initialization 
     338                           CALL dta_dyn_init( Nbb, Nnn, Naa )        ! Initialization for the dynamics 
     339#endif 
     340 
     341                           CALL     trc_init( Nbb, Nnn, Naa )        ! Passive tracers initialization 
    325342                           CALL dia_ptr_init   ! Poleward TRansports initialization 
    326343                            
     
    510527   END SUBROUTINE nemo_set_cfctl 
    511528 
    512    SUBROUTINE istate_init 
     529   SUBROUTINE istate_init( Kmm, Kaa ) 
    513530      !!---------------------------------------------------------------------- 
    514531      !!                   ***  ROUTINE istate_init  *** 
     
    516533      !! ** Purpose :   Initialization to zero of the dynamics and tracers. 
    517534      !!---------------------------------------------------------------------- 
     535      INTEGER, INTENT(in) ::   Kmm, Kaa  ! ocean time level indices 
    518536      ! 
    519537      !     now fields         !     after fields      ! 
    520       un   (:,:,:)   = 0._wp   ;   ua(:,:,:) = 0._wp   ! 
    521       vn   (:,:,:)   = 0._wp   ;   va(:,:,:) = 0._wp   ! 
    522       wn   (:,:,:)   = 0._wp   !                       ! 
    523       hdivn(:,:,:)   = 0._wp   !                       ! 
    524       tsn  (:,:,:,:) = 0._wp   !                       ! 
     538      uu   (:,:,:,Kmm)   = 0._wp   ;   uu(:,:,:,Kaa) = 0._wp   ! 
     539      vv   (:,:,:,Kmm)   = 0._wp   ;   vv(:,:,:,Kaa) = 0._wp   ! 
     540      ww   (:,:,:)   = 0._wp   !                       ! 
     541      hdiv (:,:,:)   = 0._wp   !                       ! 
     542      ts  (:,:,:,:,Kmm) = 0._wp   !                       ! 
    525543      ! 
    526544      rhd  (:,:,:) = 0.e0 
Note: See TracChangeset for help on using the changeset viewer.