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 11001 – NEMO

Changeset 11001


Ignore:
Timestamp:
2019-05-20T14:19:17+02:00 (5 years ago)
Author:
davestorkey
Message:

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : C1D and step.F90. Passes SETTE. Compiles with key_c1d.

Location:
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/C1D/dtauvd.F90

    r10068 r11001  
    117117 
    118118 
    119    SUBROUTINE dta_uvd( kt, puvd ) 
     119   SUBROUTINE dta_uvd( kt, Kmm, puvd ) 
    120120      !!---------------------------------------------------------------------- 
    121121      !!                   ***  ROUTINE dta_uvd  *** 
     
    133133      !!---------------------------------------------------------------------- 
    134134      INTEGER                           , INTENT(in   ) ::   kt     ! ocean time-step 
     135      INTEGER                           , INTENT(in   ) ::   Kmm    ! time level index 
    135136      REAL(wp), DIMENSION(jpi,jpj,jpk,2), INTENT(  out) ::   puvd   ! U & V current data 
    136137      ! 
     
    160161            DO ji = 1, jpi                ! determines the interpolated U & V current profiles at each (i,j) point 
    161162               DO jk = 1, jpk 
    162                   zl = gdept_n(ji,jj,jk) 
     163                  zl = gdept(ji,jj,jk,Kmm) 
    163164                  IF    ( zl < gdept_1d(1  ) ) THEN          ! extrapolate above the first level of data 
    164165                     zup(jk) =  puvd(ji,jj,1    ,1) 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/C1D/dyncor_c1d.F90

    r10068 r11001  
    5656 
    5757 
    58    SUBROUTINE dyn_cor_c1d( kt ) 
     58   SUBROUTINE dyn_cor_c1d( kt, Kmm, puu, pvv, Krhs ) 
    5959      !!---------------------------------------------------------------------- 
    6060      !!                   ***  ROUTINE dyn_cor_c1d  *** 
     
    6363      !!               the general trend of the momentum equation in 1D case. 
    6464      !!---------------------------------------------------------------------- 
    65       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
     65      INTEGER                             , INTENT(in   ) ::   kt        ! ocean time-step index 
     66      INTEGER                             , INTENT(in   ) ::   Kmm, Krhs ! ocean time level indices 
     67      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv  ! ocean velocities and RHS of momentum equation 
    6668      !! 
    6769      INTEGER ::   ji, jj, jk   ! dummy loop indices 
     
    7880            DO jj = 2, jpjm1 
    7981               DO ji = fs_2, fs_jpim1   ! vector opt. 
    80                   ua(ji,jj,jk) = ua(ji,jj,jk) + ff_t(ji,jj) * (vn(ji,jj,jk) + vsd(ji,jj,jk)) 
    81                   va(ji,jj,jk) = va(ji,jj,jk) - ff_t(ji,jj) * (un(ji,jj,jk) + usd(ji,jj,jk)) 
     82                  puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ff_t(ji,jj) * (pvv(ji,jj,jk,Kmm) + vsd(ji,jj,jk)) 
     83                  pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ff_t(ji,jj) * (puu(ji,jj,jk,Kmm) + usd(ji,jj,jk)) 
    8284               END DO 
    8385            END DO 
     
    8789            DO jj = 2, jpjm1 
    8890               DO ji = fs_2, fs_jpim1   ! vector opt. 
    89                   ua(ji,jj,jk) = ua(ji,jj,jk) + ff_t(ji,jj) * vn(ji,jj,jk) 
    90                   va(ji,jj,jk) = va(ji,jj,jk) - ff_t(ji,jj) * un(ji,jj,jk) 
     91                  puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ff_t(ji,jj) * pvv(ji,jj,jk,Kmm) 
     92                  pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ff_t(ji,jj) * puu(ji,jj,jk,Kmm) 
    9193               END DO 
    9294            END DO 
     
    9597       
    9698      ! 
    97       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' cor  - Ua: ', mask1=umask,  & 
    98          &                       tab3d_2=va, clinfo2=' Va: '       , mask2=vmask ) 
     99      IF(ln_ctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' cor  - Ua: ', mask1=umask,  & 
     100         &                       tab3d_2=pvv(:,:,:,Krhs), clinfo2=' Va: '       , mask2=vmask ) 
    99101      ! 
    100102   END SUBROUTINE dyn_cor_c1d 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/C1D/dyndmp.F90

    r10425 r11001  
    130130 
    131131 
    132    SUBROUTINE dyn_dmp( kt ) 
     132   SUBROUTINE dyn_dmp( kt, Kbb, Kmm, puu, pvv, Krhs ) 
    133133      !!---------------------------------------------------------------------- 
    134134      !!                   ***  ROUTINE dyn_dmp  *** 
     
    140140      !! ** Method  :   Compute Newtonian damping towards u_dta and v_dta  
    141141      !!      and add to the general momentum trends: 
    142       !!                     ua = ua + resto_uv * (u_dta - ub) 
    143       !!                     va = va + resto_uv * (v_dta - vb) 
     142      !!                     puu(Krhs) = puu(Krhs) + resto_uv * (u_dta - puu(Kbb)) 
     143      !!                     pvv(Krhs) = pvv(Krhs) + resto_uv * (v_dta - pvv(Kbb)) 
    144144      !!      The trend is computed either throughout the water column 
    145145      !!      (nn_zdmp=0), where the vertical mixing is weak (nn_zdmp=1) or 
    146146      !!      below the well mixed layer (nn_zdmp=2) 
    147147      !! 
    148       !! ** Action  : - (ua,va)   momentum trends updated with the damping trend 
    149       !!---------------------------------------------------------------------- 
    150       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     148      !! ** Action  : - (puu(:,:,:,Krhs),pvv(:,:,:,Krhs))   momentum trends updated with the damping trend 
     149      !!---------------------------------------------------------------------- 
     150      INTEGER                             , INTENT(in   ) ::   kt             ! ocean time-step index 
     151      INTEGER                             , INTENT(in   ) ::   Kbb, Kmm, Krhs ! ocean time level indices 
     152      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv       ! ocean velocities and RHS of momentum equation 
    151153      !! 
    152154      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     
    159161      ! 
    160162      !                           !==   read and interpolate U & V current data at kt   ==! 
    161       CALL dta_uvd( kt, zuv_dta ) !!! NOTE: This subroutine must be altered for use outside 
     163      CALL dta_uvd( kt, Kmm, zuv_dta ) !!! NOTE: This subroutine must be altered for use outside 
    162164                                  !!!       the C1D context (use of U,V grid variables) 
    163165      ! 
     
    168170            DO jj = 2, jpjm1 
    169171               DO ji = fs_2, fs_jpim1   ! vector opt. 
    170                   zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - ub(ji,jj,jk) ) 
    171                   zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - vb(ji,jj,jk) ) 
    172                   ua(ji,jj,jk) = ua(ji,jj,jk) + zua 
    173                   va(ji,jj,jk) = va(ji,jj,jk) + zva 
     172                  zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - puu(ji,jj,jk,Kbb) ) 
     173                  zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - pvv(ji,jj,jk,Kbb) ) 
     174                  puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zua 
     175                  pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zva 
    174176                  utrdmp(ji,jj,jk) = zua           ! save the trends 
    175177                  vtrdmp(ji,jj,jk) = zva       
     
    183185               DO ji = fs_2, fs_jpim1   ! vector opt. 
    184186                  IF( avt(ji,jj,jk) <= avt_c ) THEN 
    185                      zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - ub(ji,jj,jk) ) 
    186                      zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - vb(ji,jj,jk) ) 
     187                     zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - puu(ji,jj,jk,Kbb) ) 
     188                     zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - pvv(ji,jj,jk,Kbb) ) 
    187189                  ELSE 
    188190                     zua = 0._wp 
    189191                     zva = 0._wp   
    190192                  ENDIF 
    191                   ua(ji,jj,jk) = ua(ji,jj,jk) + zua 
    192                   va(ji,jj,jk) = va(ji,jj,jk) + zva 
     193                  puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zua 
     194                  pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zva 
    193195                  utrdmp(ji,jj,jk) = zua           ! save the trends 
    194196                  vtrdmp(ji,jj,jk) = zva 
     
    201203            DO jj = 2, jpjm1 
    202204               DO ji = fs_2, fs_jpim1   ! vector opt. 
    203                   IF( gdept_n(ji,jj,jk) >= hmlp (ji,jj) ) THEN 
    204                      zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - ub(ji,jj,jk) ) 
    205                      zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - vb(ji,jj,jk) ) 
     205                  IF( gdept(ji,jj,jk,Kmm) >= hmlp (ji,jj) ) THEN 
     206                     zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - puu(ji,jj,jk,Kbb) ) 
     207                     zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - pvv(ji,jj,jk,Kbb) ) 
    206208                  ELSE 
    207209                     zua = 0._wp 
    208210                     zva = 0._wp   
    209211                  ENDIF 
    210                   ua(ji,jj,jk) = ua(ji,jj,jk) + zua 
    211                   va(ji,jj,jk) = va(ji,jj,jk) + zva 
     212                  puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zua 
     213                  pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zva 
    212214                  utrdmp(ji,jj,jk) = zua           ! save the trends 
    213215                  vtrdmp(ji,jj,jk) = zva 
     
    219221      ! 
    220222      !                           ! Control print 
    221       IF( ln_ctl   )   CALL prt_ctl( tab3d_1=ua(:,:,:), clinfo1=' dmp  - Ua: ', mask1=umask,   & 
    222          &                           tab3d_2=va(:,:,:), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     223      IF( ln_ctl   )   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' dmp  - Ua: ', mask1=umask,   & 
     224         &                           tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    223225      ! 
    224226      ! 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/C1D/step_c1d.F90

    r10998 r11001  
    6969 
    7070      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    71       ! Ocean physics update                (ua, va, ta, sa used as workspace) 
     71      ! Ocean physics update         
    7272      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    73                          CALL eos_rab( tsb, rab_b )   ! before local thermal/haline expension ratio at T-points 
    74                          CALL eos_rab( tsn, rab_n )   ! now    local thermal/haline expension ratio at T-points 
    75                          CALL bn2( tsb, rab_b, rn2b ) ! before Brunt-Vaisala frequency 
    76                          CALL bn2( tsn, rab_n, rn2 ) ! now    Brunt-Vaisala frequency 
     73                         CALL eos_rab( ts(:,:,:,:,Nbb), rab_b, Nnn )  ! before local thermal/haline expension ratio at T-points 
     74                         CALL eos_rab( ts(:,:,:,:,Nnn), rab_n, Nnn )  ! now    local thermal/haline expension ratio at T-points 
     75                         CALL bn2( ts(:,:,:,:,Nbb), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency 
     76                         CALL bn2( ts(:,:,:,:,Nnn), rab_n, rn2 , Nnn ) ! now    Brunt-Vaisala frequency 
    7777       
    7878      !  VERTICAL PHYSICS 
    79                          CALL zdf_phy( kstp )         ! vertical physics update (bfr, avt, avs, avm + MLD) 
     79                         CALL zdf_phy( kstp, Nbb, Nnn, Nrhs  )    ! vertical physics update (bfr, avt, avs, avm + MLD) 
    8080 
    8181      IF(.NOT.ln_linssh )   CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh, Naa )  ! after ssh (includes call to div_hor) 
     
    8484      IF(.NOT.ln_linssh )   CALL wzv           ( kstp, Nbb, Nnn, ww,  Naa )  ! now cross-level velocity  
    8585      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    86       ! diagnostics and outputs             (ua, va, ta, sa used as workspace) 
     86      ! diagnostics and outputs        
    8787      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    8888                         CALL dia_wri( kstp, Nnn )  ! ocean model: outputs 
     
    9494      ! Passive Tracer Model 
    9595      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    96                         CALL trc_stp( kstp )       ! time-stepping 
     96                        CALL trc_stp( kstp, Nbb, Nnn, Nrhs, Naa  )   ! time-stepping 
    9797#endif 
    9898 
    9999      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    100       ! Active tracers                              (ua, va used as workspace) 
     100      ! Active tracers                              (uu(:,:,:,Nrhs), vv(:,:,:,Nrhs) used as workspace) 
    101101      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    102                         tsa(:,:,:,:) = 0._wp       ! set tracer trends to zero 
     102                        ts(:,:,:,:,Nrhs) = 0._wp       ! set tracer trends to zero 
    103103 
    104                         CALL tra_sbc( kstp )       ! surface boundary condition 
    105       IF( ln_traqsr )   CALL tra_qsr( kstp )       ! penetrative solar radiation qsr 
    106       IF( ln_tradmp )   CALL tra_dmp( kstp )       ! internal damping trends- tracers 
    107       IF(.NOT.ln_linssh)CALL tra_adv( kstp )       ! horizontal & vertical advection 
    108       IF( ln_zdfosm  )  CALL tra_osm( kstp )       ! OSMOSIS non-local tracer fluxes 
    109                         CALL tra_zdf( kstp )      ! vertical mixing 
    110                         CALL eos( tsn, rhd, rhop, gdept_0(:,:,:) )   ! now potential density for zdfmxl 
    111       IF( ln_zdfnpc )   CALL tra_npc( kstp )      ! applied non penetrative convective adjustment on (t,s) 
    112                         CALL tra_nxt( kstp )      ! tracer fields at next time step 
     104                        CALL tra_sbc( kstp,      Nnn, ts, Nrhs  )  ! surface boundary condition 
     105      IF( ln_traqsr )   CALL tra_qsr( kstp,      Nnn, ts, Nrhs  )  ! penetrative solar radiation qsr 
     106      IF( ln_tradmp )   CALL tra_dmp( kstp, Nbb, Nnn, ts, Nrhs  )  ! internal damping trends- tracers 
     107      IF(.NOT.ln_linssh)CALL tra_adv( kstp, Nbb, Nnn, ts, Nrhs  )  ! horizontal & vertical advection 
     108      IF( ln_zdfosm  )  CALL tra_osm( kstp, Nnn     , ts, Nrhs  )  ! OSMOSIS non-local tracer fluxes 
     109                        CALL tra_zdf( kstp, Nbb, Nnn, Nrhs, ts, Naa   ) ! vertical mixing 
     110                        CALL eos( ts(:,:,:,:,Nnn), rhd, rhop, gdept_0(:,:,:) )  ! now potential density for zdfmxl 
     111      IF( ln_zdfnpc )   CALL tra_npc( kstp,      Nnn, Nrhs, ts, Naa   ) ! applied non penetrative convective adjustment on (t,s) 
     112                        CALL tra_nxt( kstp, Nbb, Nnn, Nrhs,     Naa   ) ! tracer fields at next time step 
    113113 
    114114 
     
    117117      ! Dynamics                                    (ta, sa used as workspace) 
    118118      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    119                         ua(:,:,:) = 0._wp          ! set dynamics trends to zero 
    120                         va(:,:,:) = 0._wp 
     119                        uu(:,:,:,Nrhs) = 0._wp          ! set dynamics trends to zero 
     120                        vv(:,:,:,Nrhs) = 0._wp 
    121121 
    122       IF( ln_dyndmp )   CALL dyn_dmp    ( kstp )   ! internal damping trends- momentum 
    123                         CALL dyn_cor_c1d( kstp )   ! vorticity term including Coriolis 
    124       IF( ln_zdfosm  )  CALL dyn_osm    ( kstp )   ! OSMOSIS non-local velocity fluxes 
    125                         CALL dyn_zdf    ( kstp )   ! vertical diffusion 
    126                         CALL dyn_nxt    ( kstp )   ! lateral velocity at next time step 
    127       IF(.NOT.ln_linssh)CALL ssh_swp    ( kstp, Nbb, Nnn, Naa )   ! swap of sea surface height 
     122      IF( ln_dyndmp )   CALL dyn_dmp    ( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! internal damping trends- momentum 
     123                        CALL dyn_cor_c1d( kstp,      Nnn      , uu, vv, Nrhs )  ! vorticity term including Coriolis 
     124      IF( ln_zdfosm  )  CALL dyn_osm    ( kstp,      Nnn      , uu, vv, Nrhs )  ! OSMOSIS non-local velocity fluxes 
     125                        CALL dyn_zdf    ( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa  )  ! vertical diffusion 
     126                        CALL dyn_nxt    ( kstp, Nbb, Nnn, Naa )                 ! lateral velocity at next time step 
     127      IF(.NOT.ln_linssh)CALL ssh_swp    ( kstp, Nbb, Nnn, Naa )                 ! swap of sea surface height 
    128128 
    129       IF(.NOT.ln_linssh)CALL dom_vvl_sf_swp( kstp, Nbb, Nnn, Naa )! swap of vertical scale factors 
     129      IF(.NOT.ln_linssh)CALL dom_vvl_sf_swp( kstp, Nbb, Nnn, Naa )              ! swap of vertical scale factors 
    130130 
    131131      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    132132      ! Control and restarts 
    133133      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    134                              CALL stp_ctl( kstp, indic ) 
    135       IF( kstp == nit000 )   CALL iom_close( numror )      ! close input  ocean restart file 
    136       IF( lrst_oce       )   CALL rst_write( kstp )        ! write output ocean restart file 
     134                             CALL stp_ctl( kstp, Nnn, indic ) 
     135      IF( kstp == nit000 )   CALL iom_close( numror )          ! close input  ocean restart file 
     136      IF( lrst_oce       )   CALL rst_write( kstp, Nbb, Nnn )  ! write output ocean restart file 
    137137      ! 
    138138#if defined key_iomput 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90

    r10998 r11001  
    125125      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    126126      IF( ln_sto_eos ) CALL sto_par( kstp )          ! Stochastic parameters 
    127       IF( ln_sto_eos ) CALL sto_pts( tsn  )          ! Random T/S fluctuations 
     127      IF( ln_sto_eos ) CALL sto_pts( ts(:,:,:,:,Nnn)  )          ! Random T/S fluctuations 
    128128 
    129129      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    131131      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    132132      !  THERMODYNAMICS 
    133                          CALL eos_rab( tsb, rab_b, Nnn )       ! before local thermal/haline expension ratio at T-points 
    134                          CALL eos_rab( tsn, rab_n, Nnn )       ! now    local thermal/haline expension ratio at T-points 
    135                          CALL bn2    ( tsb, rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency 
    136                          CALL bn2    ( tsn, rab_n, rn2, Nnn  ) ! now    Brunt-Vaisala frequency 
     133                         CALL eos_rab( ts(:,:,:,:,Nbb), rab_b, Nnn )       ! before local thermal/haline expension ratio at T-points 
     134                         CALL eos_rab( ts(:,:,:,:,Nnn), rab_n, Nnn )       ! now    local thermal/haline expension ratio at T-points 
     135                         CALL bn2    ( ts(:,:,:,:,Nbb), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency 
     136                         CALL bn2    ( ts(:,:,:,:,Nnn), rab_n, rn2, Nnn  ) ! now    Brunt-Vaisala frequency 
    137137 
    138138      !  VERTICAL PHYSICS 
     
    142142      ! 
    143143      IF( l_ldfslp ) THEN                             ! slope of lateral mixing 
    144                          CALL eos( tsb, rhd, gdept_0(:,:,:) )               ! before in situ density 
     144                         CALL eos( ts(:,:,:,:,Nbb), rhd, gdept_0(:,:,:) )               ! before in situ density 
    145145 
    146146         IF( ln_zps .AND. .NOT. ln_isfcav)                                    & 
    147             &            CALL zps_hde    ( kstp, Nnn, jpts, tsb, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient 
     147            &            CALL zps_hde    ( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv,  &  ! Partial steps: before horizontal gradient 
    148148            &                                          rhd, gru , grv    )       ! of t, s, rd at the last ocean level 
    149149 
    150150         IF( ln_zps .AND.       ln_isfcav)                                                & 
    151             &            CALL zps_hde_isf( kstp, Nnn, jpts, tsb, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF) 
     151            &            CALL zps_hde_isf( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF) 
    152152            &                                          rhd, gru , grv , grui, grvi   )       ! of t, s, rd at the first ocean level 
    153153         IF( ln_traldf_triad ) THEN  
     
    169169                            CALL wzv           ( kstp, Nbb, Nnn, ww,  Naa )    ! now cross-level velocity  
    170170      IF( ln_zad_Aimp )     CALL wAimp         ( kstp,      Nnn           )  ! Adaptive-implicit vertical advection partitioning 
    171                             CALL eos    ( tsn, rhd, rhop, gdept_n(:,:,:) )  ! now in situ density for hpg computation 
     171                            CALL eos    ( ts(:,:,:,:,Nnn), rhd, rhop, gdept(:,:,:,Nnn) )  ! now in situ density for hpg computation 
    172172                             
    173173!!jc: fs simplification 
     
    176176!!                                         with previous versions using split-explicit free surface           
    177177            IF( ln_zps .AND. .NOT. ln_isfcav )                                    & 
    178                &            CALL zps_hde    ( kstp, Nnn, jpts, tsn, gtsu, gtsv,   &  ! Partial steps: before horizontal gradient 
     178               &            CALL zps_hde    ( kstp, Nnn, jpts, ts(:,:,:,:,Nnn), gtsu, gtsv,   &  ! Partial steps: before horizontal gradient 
    179179               &                                          rhd, gru , grv     )       ! of t, s, rd at the last ocean level 
    180180            IF( ln_zps .AND.       ln_isfcav )                                               & 
    181                &            CALL zps_hde_isf( kstp, Nnn, jpts, tsn, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF) 
     181               &            CALL zps_hde_isf( kstp, Nnn, jpts, ts(:,:,:,:,Nnn), gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF) 
    182182               &                                          rhd, gru , grv , grui, grvi   )       ! of t, s, rd at the first ocean level 
    183183!!jc: fs simplification 
    184184                             
    185                          ua(:,:,:) = 0._wp            ! set dynamics trends to zero 
    186                          va(:,:,:) = 0._wp 
     185                         uu(:,:,:,Nrhs) = 0._wp            ! set dynamics trends to zero 
     186                         vv(:,:,:,Nrhs) = 0._wp 
    187187 
    188188      IF(  lk_asminc .AND. ln_asmiau .AND. ln_dyninc )   & 
     
    200200                         CALL dyn_spg( kstp, Nbb, Nnn, Nrhs, uu, vv, ssh, uu_b, vv_b, Naa )  ! surface pressure gradient 
    201201 
    202                                                       ! With split-explicit free surface, since now transports have been updated and ssha as well 
     202                                                      ! With split-explicit free surface, since now transports have been updated and ssh(:,:,Nrhs) as well 
    203203      IF( ln_dynspg_ts ) THEN                         ! vertical scale factors and vertical velocity need to be updated 
    204204                            CALL div_hor       ( kstp, Nbb, Nnn )    ! Horizontal divergence  (2nd call in time-split case) 
     
    238238      ! Active tracers                               
    239239      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    240                          tsa(:,:,:,:) = 0._wp         ! set tracer trends to zero 
     240                         ts(:,:,:,:,Nrhs) = 0._wp         ! set tracer trends to zero 
    241241 
    242242      IF(  lk_asminc .AND. ln_asmiau .AND. & 
    243          & ln_trainc )   CALL tra_asm_inc   ( kstp, Nbb, Nnn, ts, Nrhs )  ! apply tracer assimilation increment 
    244                          CALL tra_sbc       ( kstp,      Nnn, ts, Nrhs )  ! surface boundary condition 
    245       IF( ln_traqsr  )   CALL tra_qsr       ( kstp,      Nnn, ts, Nrhs )  ! penetrative solar radiation qsr 
    246       IF( ln_trabbc  )   CALL tra_bbc       ( kstp,      Nnn, ts, Nrhs )  ! bottom heat flux 
    247       IF( ln_trabbl  )   CALL tra_bbl       ( kstp, Nbb, Nnn, ts, Nrhs )  ! advective (and/or diffusive) bottom boundary layer scheme 
    248       IF( ln_tradmp  )   CALL tra_dmp       ( kstp, Nbb, Nnn, ts, Nrhs )  ! internal damping trends 
    249       IF( ln_bdy     )   CALL bdy_tra_dmp   ( kstp, Nbb,      ts, Nrhs )  ! bdy damping trends 
     243         & ln_trainc )   CALL tra_asm_inc( kstp, Nbb, Nnn, ts, Nrhs )  ! apply tracer assimilation increment 
     244                         CALL tra_sbc    ( kstp,      Nnn, ts, Nrhs )  ! surface boundary condition 
     245      IF( ln_traqsr  )   CALL tra_qsr    ( kstp,      Nnn, ts, Nrhs )  ! penetrative solar radiation qsr 
     246      IF( ln_trabbc  )   CALL tra_bbc    ( kstp,      Nnn, ts, Nrhs )  ! bottom heat flux 
     247      IF( ln_trabbl  )   CALL tra_bbl    ( kstp, Nbb, Nnn, ts, Nrhs )  ! advective (and/or diffusive) bottom boundary layer scheme 
     248      IF( ln_tradmp  )   CALL tra_dmp    ( kstp, Nbb, Nnn, ts, Nrhs )  ! internal damping trends 
     249      IF( ln_bdy     )   CALL bdy_tra_dmp( kstp, Nbb,      ts, Nrhs )  ! bdy damping trends 
    250250#if defined key_agrif 
    251251      IF(.NOT. Agrif_Root())  &  
    252252               &         CALL Agrif_Sponge_tra        ! tracers sponge 
    253253#endif 
    254                          CALL tra_adv( kstp, Nbb, Nnn      , ts, Nrhs )  ! hor. + vert. advection  ==> RHS 
    255       IF( ln_zdfosm  )   CALL tra_osm( kstp,      Nnn      , ts, Nrhs )  ! OSMOSIS non-local tracer fluxes ==> RHS 
     254                         CALL tra_adv    ( kstp, Nbb, Nnn, ts, Nrhs )  ! hor. + vert. advection ==> RHS 
     255      IF( ln_zdfosm  )   CALL tra_osm    ( kstp, Nnn     , ts, Nrhs )  ! OSMOSIS non-local tracer fluxes ==> RHS 
    256256      IF( lrst_oce .AND. ln_zdfosm ) & 
    257            &             CALL osm_rst( kstp,      Nnn, 'WRITE' )         ! write OSMOSIS outputs + wn (so must do here) to restarts 
    258                          CALL tra_ldf( kstp, Nbb, Nnn      , ts, Nrhs )  ! lateral mixing 
     257           &             CALL osm_rst    ( kstp, Nnn, 'WRITE' )        ! write OSMOSIS outputs + ww (so must do here) to restarts 
     258                         CALL tra_ldf    ( kstp, Nbb, Nnn, ts, Nrhs )  ! lateral mixing 
    259259 
    260260!!gm : why CALL to dia_ptr has been moved here??? (use trends info?) 
    261261      IF( ln_diaptr  )   CALL dia_ptr( Nnn )                 ! Poleward adv/ldf TRansports diagnostics 
    262262!!gm 
    263                          CALL tra_zdf( kstp, Nbb, Nnn, Nrhs, ts, Naa  )  ! vert. mixing & after tracer   ==> after 
    264       IF( ln_zdfnpc  )   CALL tra_npc( kstp,      Nnn, Nrhs, ts, Naa  )  ! update after fields by non-penetrative convection 
     263                         CALL tra_zdf    ( kstp, Nbb, Nnn, Nrhs, ts, Naa  )  ! vert. mixing & after tracer  ==> after 
     264      IF( ln_zdfnpc  )   CALL tra_npc    ( kstp,      Nnn, Nrhs, ts, Naa  )  ! update after fields by non-penetrative convection 
    265265 
    266266      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    280280!!    place. 
    281281!!  
    282 !!jc2: dynnxt must be the latest call. e3t_b are indeed updated in that routine 
     282!!jc2: dynnxt must be the latest call. e3t(:,:,:,Nbb) are indeed updated in that routine 
    283283                         CALL tra_nxt       ( kstp, Nbb, Nnn, Nrhs, Naa )  ! finalize (bcs) tracer fields at next time step and swap 
    284284                         CALL dyn_nxt       ( kstp, Nbb, Nnn, Naa  )  ! finalize (bcs) velocities at next time step and swap (always called after tra_nxt) 
Note: See TracChangeset for help on using the changeset viewer.