Changeset 10806


Ignore:
Timestamp:
2019-03-27T17:55:22+01:00 (18 months ago)
Author:
davestorkey
Message:

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps branch: Latest updates. Make sure all time-dependent 3D variables are converted in previously modified modules.

Location:
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps
Files:
24 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/cfgs/SHARED/field_def_nemo-oce.xml

    r10645 r10806  
    6262        <field id="mldr10_1max"  long_name="Max of Mixed Layer Depth (dsigma = 0.01 wrt 10m)"   field_ref="mldr10_1"   operation="maximum"                                                                          /> 
    6363        <field id="mldr10_1min"  long_name="Min of Mixed Layer Depth (dsigma = 0.01 wrt 10m)"   field_ref="mldr10_1"   operation="minimum"                                                                          /> 
     64         <field id="mldzint_1"    long_name="Mixed Layer Depth interpolated"                     standard_name="ocean_mixed_layer_thickness"                                                       unit="m"          /> 
     65         <field id="mldzint_2"    long_name="Mixed Layer Depth interpolated"                     standard_name="ocean_mixed_layer_thickness"                                                       unit="m"          /> 
     66         <field id="mldzint_3"    long_name="Mixed Layer Depth interpolated"                     standard_name="ocean_mixed_layer_thickness"                                                       unit="m"          /> 
     67         <field id="mldzint_4"    long_name="Mixed Layer Depth interpolated"                     standard_name="ocean_mixed_layer_thickness"                                                       unit="m"          /> 
     68         <field id="mldzint_5"    long_name="Mixed Layer Depth interpolated"                     standard_name="ocean_mixed_layer_thickness"                                                       unit="m"          /> 
     69         <field id="mldhtc_1"     long_name="Mixed Layer Depth integrated heat content"          standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       /> 
     70         <field id="mldhtc_2"     long_name="Mixed Layer Depth integrated heat content"          standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       /> 
     71         <field id="mldhtc_3"     long_name="Mixed Layer Depth integrated heat content"          standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       /> 
     72         <field id="mldhtc_4"     long_name="Mixed Layer Depth integrated heat content"          standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       /> 
     73         <field id="mldhtc_5"     long_name="Mixed Layer Depth integrated heat content"          standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       /> 
    6474        <field id="heatc"        long_name="Heat content vertically integrated"                 standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       /> 
    6575        <field id="saltc"        long_name="Salt content vertically integrated"                                                                                                                   unit="1e-3*kg/m2" /> 
     
    7787        <!-- variables available with MLE --> 
    7888        <field id="Lf_NHpf"      long_name="MLE: Lf = N H / f"   unit="m" /> 
    79  
    80         <!-- next variables available with ln_zad_Aimp=.true. --> 
    81         <field id="Courant"      long_name="Courant number"                                                                 unit="#"   grid_ref="grid_T_3D" /> 
    82         <field id="wimp"         long_name="Implicit vertical velocity"                                                     unit="m/s" grid_ref="grid_T_3D" /> 
    83         <field id="wexp"         long_name="Explicit vertical velocity"                                                     unit="m/s" grid_ref="grid_T_3D" /> 
    84         <field id="wi_cff"       long_name="Fraction of implicit vertical velocity"                                         unit="#"   grid_ref="grid_T_3D" /> 
    8589 
    8690        <!-- next variables available with key_diahth --> 
     
    351355        <field id="uoce"         long_name="ocean current along i-axis"                             standard_name="sea_water_x_velocity"        unit="m/s"        grid_ref="grid_U_3D" /> 
    352356        <field id="uoce_e3u"     long_name="ocean current along i-axis  (thickness weighted)"                                                   unit="m/s"        grid_ref="grid_U_3D"  > uoce * e3u </field> 
     357        <field id="uoce2_e3u"    long_name="ocean current along i-axis squared (thickness weighted)"                                            unit="m3/s2"      grid_ref="grid_U_3D"  > uoce * uoce * e3u </field> 
    353358        <field id="ssu"          long_name="ocean surface current along i-axis"                                                                 unit="m/s"                             /> 
    354359        <field id="sbu"          long_name="ocean bottom current along i-axis"                                                                  unit="m/s"                             /> 
     
    405410        <field id="voce"         long_name="ocean current along j-axis"                             standard_name="sea_water_y_velocity"        unit="m/s"        grid_ref="grid_V_3D" /> 
    406411        <field id="voce_e3v"     long_name="ocean current along j-axis  (thickness weighted)"                                                   unit="m/s"        grid_ref="grid_V_3D"  > voce * e3v </field> 
     412        <field id="voce2_e3v"    long_name="ocean current along j-axis squared (thickness weighted)"                                            unit="m3/s2"      grid_ref="grid_V_3D"  > voce * voce * e3v </field> 
    407413        <field id="ssv"          long_name="ocean surface current along j-axis"                                                                 unit="m/s"                             /> 
    408414        <field id="sbv"          long_name="ocean bottom current along j-axis"                                                                  unit="m/s"                             /> 
     
    975981 
    976982    <field_group id="25h_grid_W" grid_ref="grid_W_3D" operation="instant"> 
    977       <field id="vovecrtz25h"         name="k current 25h mean"                 unit="m/s"      /> 
     983      <field id="vomecrtz25h"         name="k current 25h mean"                 unit="m/s"      /> 
    978984      <field id="avt25h"              name="vertical diffusivity25h mean"       unit="m2/s" /> 
    979985      <field id="avm25h"              name="vertical viscosity 25h mean"        unit="m2/s" /> 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/dom_oce.F90

    r10787 r10806  
    154154   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:,:)   ::  gde3w   
    155155    
    156    !                                                      !  reference depths of water column 
    157    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht_0  !: t-depth              [m] 
    158    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_0  !: u-depth              [m] 
    159    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hv_0  !: v-depth              [m] 
    160    !                                                      !  time-dependent depths of water column 
    161    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:,:) ::  hu, hv, r1_hu, r1_hv   
    162    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:)   ::  ht 
     156   !                                                      !  ref. ! before  !   now   !  after  ! 
     157   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht_0            ,    ht_n             !: t-depth              [m] 
     158   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_0  ,    hu_b ,    hu_n ,    hu_a   !: u-depth              [m] 
     159   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hv_0  ,    hv_b ,    hv_n ,    hv_a   !: v-depth              [m] 
     160   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::           r1_hu_b , r1_hu_n , r1_hu_a   !: inverse of u-depth [1/m] 
     161   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::           r1_hv_b , r1_hv_n , r1_hv_a   !: inverse of v-depth [1/m] 
    163162 
    164163   !! TEMPORARY POINTERS FOR DEVELOPMENT ONLY 
     
    173172   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:) ::    gdepw_b , gdepw_n           !: w- depth              [m] 
    174173   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:) ::    gde3w_n                     !: w- depth (sum of e3w) [m] 
    175    REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:)   ::    ht_n                        !: t-depth              [m] 
    176    REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:)   ::    hu_b ,    hu_n ,    hu_a    !: u-depth              [m] 
    177    REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:)   ::    hv_b ,    hv_n ,    hv_a    !: v-depth              [m] 
    178    REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:)   :: r1_hu_b , r1_hu_n , r1_hu_a    !: inverse of u-depth [1/m] 
    179    REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:)   :: r1_hv_b , r1_hv_n , r1_hv_a    !: inverse of v-depth [1/m] 
    180174   !! TEMPORARY POINTERS FOR DEVELOPMENT ONLY 
    181175 
     
    295289         &      e3f  (jpi,jpj,jpk),     STAT=ierr(5) )                        
    296290         ! 
    297       ALLOCATE( ht_0(jpi,jpj) ,    hu_0(jpi,jpj)     ,    hv_0(jpi,jpj) ,                                           & 
    298          &      ht  (jpi,jpj) ,    hu  (jpi,jpj,jpt) ,    hv  (jpi,jpj,jpt) ,                                       & 
    299          &                      r1_hu  (jpi,jpj,jpt) , r1_hv  (jpi,jpj,jpt) , STAT=ierr(6) ) 
     291      ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj) , hv_0(jpi,jpj) ,                                           & 
     292         &                      hu_b(jpi,jpj) , hv_b(jpi,jpj) , r1_hu_b(jpi,jpj) , r1_hv_b(jpi,jpj) ,     & 
     293         &      ht_n(jpi,jpj) , hu_n(jpi,jpj) , hv_n(jpi,jpj) , r1_hu_n(jpi,jpj) , r1_hv_n(jpi,jpj) ,     & 
     294         &                      hu_a(jpi,jpj) , hv_a(jpi,jpj) , r1_hu_a(jpi,jpj) , r1_hv_a(jpi,jpj) , STAT=ierr(6)  ) 
    300295         ! 
    301296         ! 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf.F90

    r10068 r10806  
    4343CONTAINS 
    4444 
    45    SUBROUTINE dyn_ldf( kt ) 
     45   SUBROUTINE dyn_ldf( kt, ktlev1, ktlev2, pu_rhs, pv_rhs ) 
    4646      !!---------------------------------------------------------------------- 
    4747      !!                  ***  ROUTINE dyn_ldf  *** 
     
    4949      !! ** Purpose :   compute the lateral ocean dynamics physics. 
    5050      !!---------------------------------------------------------------------- 
    51       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     51      INTEGER, INTENT(in) ::   kt               ! ocean time-step index 
     52      INTEGER, INTENT(in) ::   ktlev1, ktlev2   ! time level index for source terms 
     53      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pu_rhs, pv_rhs ! momentum trends 
    5254      ! 
    5355      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv 
     
    5860      IF( l_trddyn )   THEN                      ! temporary save of momentum trends 
    5961         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 
    60          ztrdu(:,:,:) = ua(:,:,:)  
    61          ztrdv(:,:,:) = va(:,:,:)  
     62         ztrdu(:,:,:) = pu_rhs(:,:,:)  
     63         ztrdv(:,:,:) = pv_rhs(:,:,:)  
    6264      ENDIF 
    6365 
    6466      SELECT CASE ( nldf_dyn )                   ! compute lateral mixing trend and add it to the general trend 
    6567      ! 
    66       CASE ( np_lap   )    ;   CALL dyn_ldf_lap( kt, ub, vb, ua, va, 1 )      ! iso-level    laplacian 
     68      CASE ( np_lap   )    ;   CALL dyn_ldf_lap( kt, ktlev1, ktlev2, uu(:,:,:,ktlev1), vv(:,:,:,ktlev1), pu_rhs, pv_rhs, 1 )      ! iso-level    laplacian 
    6769      CASE ( np_lap_i )    ;   CALL dyn_ldf_iso( kt )                         ! rotated      laplacian 
    68       CASE ( np_blp   )    ;   CALL dyn_ldf_blp( kt, ub, vb, ua, va    )      ! iso-level bi-laplacian 
     70      CASE ( np_blp   )    ;   CALL dyn_ldf_blp( kt, ktlev1, ktlev2, uu(:,:,:,ktlev1), vv(:,:,:,ktlev1), pu_rhs, pv_rhs    )      ! iso-level bi-laplacian 
    6971      ! 
    7072      END SELECT 
    7173 
    7274      IF( l_trddyn ) THEN                        ! save the horizontal diffusive trends for further diagnostics 
    73          ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    74          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     75         ztrdu(:,:,:) = pu_rhs(:,:,:) - ztrdu(:,:,:) 
     76         ztrdv(:,:,:) = pv_rhs(:,:,:) - ztrdv(:,:,:) 
    7577         CALL trd_dyn( ztrdu, ztrdv, jpdyn_ldf, kt ) 
    7678         DEALLOCATE ( ztrdu , ztrdv ) 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf_lap_blp.F90

    r10425 r10806  
    3535CONTAINS 
    3636 
    37    SUBROUTINE dyn_ldf_lap( kt, pub, pvb, pua, pva, kpass ) 
     37   SUBROUTINE dyn_ldf_lap( kt, ktlev1, ktlev2, pu, pv, pu_rhs, pva_rhs, kpass ) 
    3838      !!---------------------------------------------------------------------- 
    3939      !!                     ***  ROUTINE dyn_ldf_lap  *** 
     
    4545      !!      writen as :   grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) )  
    4646      !! 
    47       !! ** Action : - pua, pva increased by the harmonic operator applied on pub, pvb. 
     47      !! ** Action : - pu_rhs, pva_rhs increased by the harmonic operator applied on pu, pv. 
    4848      !!---------------------------------------------------------------------- 
    49       INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index 
     49      INTEGER                         , INTENT(in   ) ::   kt              ! ocean time-step index 
     50      INTEGER                         , INTENT(in   ) ::   ktlev1, ktlev2  ! time level index for scale factors 
    5051      INTEGER                         , INTENT(in   ) ::   kpass      ! =1/2 first or second passage 
    51       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity  [m/s] 
    52       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! velocity trend   [m/s2] 
     52      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu, pv   ! before velocity  [m/s] 
     53      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pva_rhs   ! velocity trend   [m/s2] 
    5354      ! 
    5455      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     
    7677!!gm open question here : e3f  at before or now ?    probably now... 
    7778!!gm note that ahmf has already been multiplied by fmask 
    78                zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f_n(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       & 
    79                   &     * (  e2v(ji  ,jj-1) * pvb(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pvb(ji-1,jj-1,jk)  & 
    80                   &        - e1u(ji-1,jj  ) * pub(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pub(ji-1,jj-1,jk)  ) 
     79               zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       & 
     80                  &     * (  e2v(ji  ,jj-1) * pv(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk)  & 
     81                  &        - e1u(ji-1,jj  ) * pu(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk)  ) 
    8182               !                                      ! ahm * div        (computed from 2 to jpi/jpj) 
    8283!!gm note that ahmt has already been multiplied by tmask 
    83                zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk)                                         & 
    84                   &     * (  e2u(ji,jj)*e3u_b(ji,jj,jk) * pub(ji,jj,jk) - e2u(ji-1,jj)*e3u_b(ji-1,jj,jk) * pub(ji-1,jj,jk)  & 
    85                   &        + e1v(ji,jj)*e3v_b(ji,jj,jk) * pvb(ji,jj,jk) - e1v(ji,jj-1)*e3v_b(ji,jj-1,jk) * pvb(ji,jj-1,jk)  ) 
     84               zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,ktlev1)                                         & 
     85                  &     * (  e2u(ji,jj)*e3u(ji,jj,jk,ktlev1) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,ktlev1) * pu(ji-1,jj,jk)  & 
     86                  &        + e1v(ji,jj)*e3v(ji,jj,jk,ktlev1) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,ktlev1) * pv(ji,jj-1,jk)  ) 
    8687            END DO   
    8788         END DO   
     
    8990         DO jj = 2, jpjm1                             ! - curl( curl) + grad( div ) 
    9091            DO ji = fs_2, fs_jpim1   ! vector opt. 
    91                pua(ji,jj,jk) = pua(ji,jj,jk) + zsign * (                                                 & 
    92                   &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk)   & 
     92               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * (                                                 & 
     93                  &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,ktlev2)   & 
    9394                  &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                     ) 
    9495                  ! 
    95                pva(ji,jj,jk) = pva(ji,jj,jk) + zsign * (                                                 & 
    96                   &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk)   & 
     96               pva_rhs(ji,jj,jk) = pva_rhs(ji,jj,jk) + zsign * (                                                 & 
     97                  &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,ktlev2)   & 
    9798                  &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                     ) 
    9899            END DO 
     
    105106 
    106107 
    107    SUBROUTINE dyn_ldf_blp( kt, pub, pvb, pua, pva ) 
     108   SUBROUTINE dyn_ldf_blp( kt, ktlev1, ktlev2, pu, pv, pu_rhs, pva_rhs ) 
    108109      !!---------------------------------------------------------------------- 
    109110      !!                 ***  ROUTINE dyn_ldf_blp  *** 
     
    116117      !!      It is computed by two successive calls to dyn_ldf_lap routine 
    117118      !! 
    118       !! ** Action :   pta   updated with the before rotated bilaplacian diffusion 
     119      !! ** Action :   pt_rhs   updated with the before rotated bilaplacian diffusion 
    119120      !!---------------------------------------------------------------------- 
    120121      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index 
    121       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity fields 
    122       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! momentum trend 
     122      INTEGER                         , INTENT(in   ) ::   ktlev1, ktlev2  ! time level index for scale factors 
     123      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu, pv   ! before velocity fields 
     124      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pva_rhs   ! momentum trend 
    123125      ! 
    124126      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulap, zvlap   ! laplacian at u- and v-point 
     
    134136      zvlap(:,:,:) = 0._wp 
    135137      ! 
    136       CALL dyn_ldf_lap( kt, pub, pvb, zulap, zvlap, 1 )   ! rotated laplacian applied to ptb (output in zlap) 
     138      CALL dyn_ldf_lap( kt, ktlev1, ktlev2, pu, pv, zulap, zvlap, 1 )   ! rotated laplacian applied to pt (output in zlap) 
    137139      ! 
    138140      CALL lbc_lnk_multi( 'dynldf_lap_blp', zulap, 'U', -1., zvlap, 'V', -1. )             ! Lateral boundary conditions 
    139141      ! 
    140       CALL dyn_ldf_lap( kt, zulap, zvlap, pua, pva, 2 )   ! rotated laplacian applied to zlap (output in pta) 
     142      CALL dyn_ldf_lap( kt, ktlev1, ktlev2, zulap, zvlap, pu_rhs, pva_rhs, 2 )   ! rotated laplacian applied to zlap (output in pt_rhs) 
    141143      ! 
    142144   END SUBROUTINE dyn_ldf_blp 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynvor.F90

    r10802 r10806  
    119119         ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 
    120120         ! 
    121          ztrdu(:,:,:) = ua(:,:,:)            !* planetary vorticity trend (including Stokes-Coriolis force) 
    122          ztrdv(:,:,:) = va(:,:,:) 
     121         ztrdu(:,:,:) = pu_rhs(:,:,:)            !* planetary vorticity trend (including Stokes-Coriolis force) 
     122         ztrdv(:,:,:) = pv_rhs(:,:,:) 
    123123         SELECT CASE( nvor_scheme ) 
    124124         CASE( np_ENS )           ;   CALL vor_ens( kt, ktlev, ncor, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )   ! enstrophy conserving scheme 
     
    133133            IF( ln_stcor )            CALL vor_een( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs )   ! add the Stokes-Coriolis trend 
    134134         END SELECT 
    135          ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    136          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     135         ztrdu(:,:,:) = pu_rhs(:,:,:) - ztrdu(:,:,:) 
     136         ztrdv(:,:,:) = pv_rhs(:,:,:) - ztrdv(:,:,:) 
    137137         CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 
    138138         ! 
    139139         IF( n_dynadv /= np_LIN_dyn ) THEN   !* relative vorticity or metric trend (only in non-linear case) 
    140             ztrdu(:,:,:) = ua(:,:,:) 
    141             ztrdv(:,:,:) = va(:,:,:) 
     140            ztrdu(:,:,:) = pu_rhs(:,:,:) 
     141            ztrdv(:,:,:) = pv_rhs(:,:,:) 
    142142            SELECT CASE( nvor_scheme ) 
    143143            CASE( np_ENT )           ;   CALL vor_enT( kt, ktlev, nrvm, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )  ! energy conserving scheme (T-pts) 
     
    147147            CASE( np_EEN )           ;   CALL vor_een( kt, ktlev, nrvm, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )  ! energy & enstrophy scheme 
    148148            END SELECT 
    149             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    150             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     149            ztrdu(:,:,:) = pu_rhs(:,:,:) - ztrdu(:,:,:) 
     150            ztrdv(:,:,:) = pv_rhs(:,:,:) - ztrdv(:,:,:) 
    151151            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt ) 
    152152         ENDIF 
     
    410410 
    411411         IF( ln_sco ) THEN 
    412             zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk) 
     412            zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
    413413            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,ktlev) * pu(:,:,jk) 
    414414            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,ktlev) * pv(:,:,jk) 
     
    518518         ! 
    519519         IF( ln_sco ) THEN                   !==  horizontal fluxes  ==! 
    520             zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk) 
     520            zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
    521521            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,ktlev) * pu(:,:,jk) 
    522522            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,ktlev) * pv(:,:,jk) 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv.F90

    r10802 r10806  
    7575CONTAINS 
    7676 
    77    SUBROUTINE tra_adv( kt, ktlev1, ktlev2, ktlev3, pts_rhs ) 
     77   SUBROUTINE tra_adv( kt, ktlev1, ktlev2, ktlev3, kt2lev, pts_rhs ) 
    7878      !!---------------------------------------------------------------------- 
    7979      !!                  ***  ROUTINE tra_adv  *** 
     
    8484      !!---------------------------------------------------------------------- 
    8585      INTEGER, INTENT(in) ::   kt                       ! ocean time-step index 
    86       INTEGER, INTENT(in) ::   ktlev1, ktlev2, ktlev3   ! time level indices for source terms 
     86      INTEGER, INTENT(in) ::   ktlev1, ktlev2, ktlev3   ! time level indices for 3-time-level source terms 
     87      INTEGER, INTENT(in) ::   kt2lev                   ! time level index for 2-time-level source terms 
    8788      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends 
    8889      ! 
     
    153154     &                         ts(:,:,:,:,ktlev1), ts(:,:,:,:,ktlev2), pts_rhs, jpts, nn_fct_h, nn_fct_v ) 
    154155      CASE ( np_MUS )                                 ! MUSCL 
    155          CALL tra_adv_mus    ( kt, nit000, ktlev2, 'TRA', r2dt, zun, zvn, zwn, ts(:,:,:,:,ktlev1),      pts_rhs, jpts        , ln_mus_ups )  
     156         CALL tra_adv_mus    ( kt, nit000, ktlev2, kt2lev, 'TRA', r2dt, zun, zvn, zwn, ts(:,:,:,:,ktlev1),      pts_rhs, jpts        , ln_mus_ups )  
    156157      CASE ( np_UBS )                                 ! UBS 
    157158         CALL tra_adv_ubs    ( kt, nit000, ktlev2, 'TRA', r2dt, zun, zvn, zwn, ts(:,:,:,:,ktlev1), ts(:,:,:,:,ktlev2), pts_rhs, jpts        , nn_ubs_v   ) 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_cen.F90

    r10802 r10806  
    4444CONTAINS 
    4545 
    46    SUBROUTINE tra_adv_cen( kt, kit000, ktlev, cdtype, pu, pv, pwn,     & 
     46   SUBROUTINE tra_adv_cen( kt, kit000, ktlev, cdtype, pu, pv, pw,     & 
    4747      &                                               pt, pt_rhs, kjpt, kn_cen_h, kn_cen_v )  
    4848      !!---------------------------------------------------------------------- 
     
    7070      INTEGER                              , INTENT(in   ) ::   kn_cen_h        ! =2/4 (2nd or 4th order scheme) 
    7171      INTEGER                              , INTENT(in   ) ::   kn_cen_v        ! =2/4 (2nd or 4th order scheme) 
    72       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pwn   ! 3 ocean velocity components 
     72      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pw   ! 3 ocean velocity components 
    7373      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt             ! now tracer fields 
    7474      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs             ! tracer trend  
     
    151151               DO jj = 2, jpjm1 
    152152                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    153                      zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( pt(ji,jj,jk,jn) + pt(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk) 
     153                     zwz(ji,jj,jk) = 0.5 * pw(ji,jj,jk) * ( pt(ji,jj,jk,jn) + pt(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk) 
    154154                  END DO 
    155155               END DO 
     
    161161               DO jj = 2, jpjm1 
    162162                  DO ji = fs_2, fs_jpim1 
    163                      zwz(ji,jj,jk) = pwn(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk) 
     163                     zwz(ji,jj,jk) = pw(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk) 
    164164                  END DO 
    165165               END DO 
     
    172172               DO jj = 1, jpj 
    173173                  DO ji = 1, jpi 
    174                      zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn)  
     174                     zwz(ji,jj, mikt(ji,jj) ) = pw(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn)  
    175175                  END DO 
    176176               END DO    
    177177            ELSE                                   ! no ice-shelf cavities (only ocean surface) 
    178                zwz(:,:,1) = pwn(:,:,1) * pt(:,:,1,jn) 
     178               zwz(:,:,1) = pw(:,:,1) * pt(:,:,1,jn) 
    179179            ENDIF 
    180180         ENDIF 
     
    194194            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pu, pt(:,:,:,jn) ) 
    195195            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pv, pt(:,:,:,jn) ) 
    196             CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, pt(:,:,:,jn) ) 
     196            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pw, pt(:,:,:,jn) ) 
    197197         END IF 
    198198         !                                 ! "Poleward" heat and salt transports  
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_fct.F90

    r10802 r10806  
    5252CONTAINS 
    5353 
    54    SUBROUTINE tra_adv_fct( kt, kit000, ktlev1, ktlev2, ktlev3, cdtype, p2dt, pu, pv, pwn,       & 
     54   SUBROUTINE tra_adv_fct( kt, kit000, ktlev1, ktlev2, ktlev3, cdtype, p2dt, pu, pv, pw,       & 
    5555      &                                                                pt_lev1, pt_lev2, pt_rhs, kjpt, kn_fct_h, kn_fct_v ) 
    5656      !!---------------------------------------------------------------------- 
     
    7777      INTEGER                              , INTENT(in   ) ::   kn_fct_v        ! order of the FCT scheme (=2 or 4) 
    7878      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step 
    79       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pwn   ! 3 ocean velocity components 
     79      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pw   ! 3 ocean velocity components 
    8080      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt_lev1, pt_lev2        ! before and now tracer fields 
    8181      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs             ! tracer trend  
     
    139139            DO jj = 1, jpj 
    140140               DO ji = 1, jpi 
    141                   zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 
    142                   zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 
     141                  zfp_wk = pw(ji,jj,jk) + ABS( pw(ji,jj,jk) ) 
     142                  zfm_wk = pw(ji,jj,jk) - ABS( pw(ji,jj,jk) ) 
    143143                  zwz(ji,jj,jk) = 0.5 * ( zfp_wk * pt_lev1(ji,jj,jk,jn) + zfm_wk * pt_lev1(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk) 
    144144               END DO 
     
    149149               DO jj = 1, jpj 
    150150                  DO ji = 1, jpi 
    151                      zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * pt_lev1(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
     151                     zwz(ji,jj, mikt(ji,jj) ) = pw(ji,jj,mikt(ji,jj)) * pt_lev1(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
    152152                  END DO 
    153153               END DO    
    154154            ELSE                             ! no cavities: only at the ocean surface 
    155                zwz(:,:,1) = pwn(:,:,1) * pt_lev1(:,:,1,jn) 
     155               zwz(:,:,1) = pw(:,:,1) * pt_lev1(:,:,1,jn) 
    156156            ENDIF 
    157157         ENDIF 
     
    258258               DO jj = 2, jpjm1 
    259259                  DO ji = fs_2, fs_jpim1 
    260                      zwz(ji,jj,jk) =  (  pwn(ji,jj,jk) * 0.5_wp * ( pt_lev2(ji,jj,jk,jn) + pt_lev2(ji,jj,jk-1,jn) )   & 
     260                     zwz(ji,jj,jk) =  (  pw(ji,jj,jk) * 0.5_wp * ( pt_lev2(ji,jj,jk,jn) + pt_lev2(ji,jj,jk-1,jn) )   & 
    261261                        &              - zwz(ji,jj,jk)  ) * wmask(ji,jj,jk) 
    262262                  END DO 
     
    269269               DO jj = 2, jpjm1 
    270270                  DO ji = fs_2, fs_jpim1 
    271                      zwz(ji,jj,jk) = ( pwn(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 
     271                     zwz(ji,jj,jk) = ( pw(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 
    272272                  END DO 
    273273               END DO 
     
    306306               CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pu, pt_lev2(:,:,:,jn) ) 
    307307               CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pv, pt_lev2(:,:,:,jn) ) 
    308                CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, pt_lev2(:,:,:,jn) ) 
     308               CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pw, pt_lev2(:,:,:,jn) ) 
    309309            ENDIF 
    310310            !                             ! heat/salt transport 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_mus.F90

    r10802 r10806  
    5454CONTAINS 
    5555 
    56    SUBROUTINE tra_adv_mus( kt, kit000, ktlev, cdtype, p2dt, pu, pv, pwn,             & 
    57       &                                                     pt, pt_rhs, kjpt, ld_msc_ups ) 
     56   SUBROUTINE tra_adv_mus( kt, kit000, ktlev, kt2lev, cdtype, p2dt, pu, pv, pw,             & 
     57      &                                                             pt, pt_rhs, kjpt, ld_msc_ups ) 
    5858      !!---------------------------------------------------------------------- 
    5959      !!                    ***  ROUTINE tra_adv_mus  *** 
     
    7575      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index 
    7676      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index 
    77       INTEGER                              , INTENT(in   ) ::   ktlev           ! time level index for source terms 
     77      INTEGER                              , INTENT(in   ) ::   ktlev           ! time level index for 3-time-level source terms 
     78      INTEGER                              , INTENT(in   ) ::   kt2lev          ! time level index for 2-time-level source terms 
    7879      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator) 
    7980      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers 
    8081      LOGICAL                              , INTENT(in   ) ::   ld_msc_ups      ! use upstream scheme within muscl 
    8182      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step 
    82       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pwn   ! 3 ocean velocity components 
     83      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pw   ! 3 ocean velocity components 
    8384      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt             ! before tracer field 
    8485      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs             ! tracer trend  
     
    240241            DO jj = 2, jpjm1       
    241242               DO ji = fs_2, fs_jpim1   ! vector opt. 
    242                   z0w = SIGN( 0.5, pwn(ji,jj,jk+1) ) 
     243                  z0w = SIGN( 0.5, pw(ji,jj,jk+1) ) 
    243244                  zalpha = 0.5 + z0w 
    244                   zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w_n(ji,jj,jk+1) 
     245                  zw  = z0w - 0.5 * pw(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w(ji,jj,jk+1,kt2lev) 
    245246                  zzwx = pt(ji,jj,jk+1,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1) 
    246247                  zzwy = pt(ji,jj,jk  ,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk  ) 
    247                   zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk) 
     248                  zwx(ji,jj,jk+1) = pw(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk) 
    248249               END DO  
    249250            END DO 
     
    253254               DO jj = 1, jpj 
    254255                  DO ji = 1, jpi 
    255                      zwx(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn) 
     256                     zwx(ji,jj, mikt(ji,jj) ) = pw(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn) 
    256257                  END DO 
    257258               END DO    
    258259            ELSE                                      ! no cavities: only at the ocean surface 
    259                zwx(:,:,1) = pwn(:,:,1) * pt(:,:,1,jn) 
     260               zwx(:,:,1) = pw(:,:,1) * pt(:,:,1,jn) 
    260261            ENDIF 
    261262         ENDIF 
     
    269270         END DO 
    270271         !                                ! send trends for diagnostic 
    271          IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, pt(:,:,:,jn) ) 
     272         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pw, pt(:,:,:,jn) ) 
    272273         ! 
    273274      END DO                     ! end of tracer loop 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_qck.F90

    r10802 r10806  
    4747CONTAINS 
    4848 
    49    SUBROUTINE tra_adv_qck ( kt, kit000, ktlev, cdtype, p2dt, pu, pv, pwn,      & 
     49   SUBROUTINE tra_adv_qck ( kt, kit000, ktlev, cdtype, p2dt, pu, pv, pw,      & 
    5050      &                                                pt_lev1, pt_lev2, pt_rhs, kjpt ) 
    5151      !!---------------------------------------------------------------------- 
     
    9090      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers 
    9191      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step 
    92       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pwn   ! 3 ocean velocity components 
     92      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pw   ! 3 ocean velocity components 
    9393      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt_lev1, pt_lev2        ! before and now tracer fields 
    9494      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs             ! tracer trend  
     
    113113 
    114114      !        ! vertical fluxes are computed with the 2nd order centered scheme 
    115       CALL tra_adv_cen2_k( kt, ktlev, cdtype, pwn,         pt_lev2, pt_rhs, kjpt ) 
     115      CALL tra_adv_cen2_k( kt, ktlev, cdtype, pw,         pt_lev2, pt_rhs, kjpt ) 
    116116      ! 
    117117   END SUBROUTINE tra_adv_qck 
     
    351351 
    352352 
    353    SUBROUTINE tra_adv_cen2_k( kt, ktlev, cdtype, pwn,           & 
     353   SUBROUTINE tra_adv_cen2_k( kt, ktlev, cdtype, pw,           & 
    354354     &                                    pt_lev2, pt_rhs, kjpt ) 
    355355      !!---------------------------------------------------------------------- 
     
    360360      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator) 
    361361      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers 
    362       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pwn      ! vertical velocity  
     362      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pw      ! vertical velocity  
    363363      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt_lev2      ! before and now tracer fields 
    364364      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs      ! tracer trend  
     
    378378            DO jj = 2, jpjm1 
    379379               DO ji = fs_2, fs_jpim1   ! vector opt. 
    380                   zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( pt_lev2(ji,jj,jk-1,jn) + pt_lev2(ji,jj,jk,jn) ) * wmask(ji,jj,jk) 
     380                  zwz(ji,jj,jk) = 0.5 * pw(ji,jj,jk) * ( pt_lev2(ji,jj,jk-1,jn) + pt_lev2(ji,jj,jk,jn) ) * wmask(ji,jj,jk) 
    381381               END DO 
    382382            END DO 
     
    386386               DO jj = 1, jpj 
    387387                  DO ji = 1, jpi 
    388                      zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * pt_lev2(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
     388                     zwz(ji,jj, mikt(ji,jj) ) = pw(ji,jj,mikt(ji,jj)) * pt_lev2(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
    389389                  END DO 
    390390               END DO    
    391391            ELSE                                   ! no ocean cavities (only ocean surface) 
    392                zwz(:,:,1) = pwn(:,:,1) * pt_lev2(:,:,1,jn) 
     392               zwz(:,:,1) = pw(:,:,1) * pt_lev2(:,:,1,jn) 
    393393            ENDIF 
    394394         ENDIF 
     
    403403         END DO 
    404404         !                                 ! Send trends for diagnostic 
    405          IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, pt_lev2(:,:,:,jn) ) 
     405         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pw, pt_lev2(:,:,:,jn) ) 
    406406         ! 
    407407      END DO 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_ubs.F90

    r10802 r10806  
    4646CONTAINS 
    4747 
    48    SUBROUTINE tra_adv_ubs( kt, kit000, ktlev, cdtype, p2dt, pu, pv, pwn,          & 
     48   SUBROUTINE tra_adv_ubs( kt, kit000, ktlev, cdtype, p2dt, pu, pv, pw,          & 
    4949      &                                                pt_lev1, pt_lev2, pt_rhs, kjpt, kn_ubs_v ) 
    5050      !!---------------------------------------------------------------------- 
     
    9191      INTEGER                              , INTENT(in   ) ::   kn_ubs_v        ! number of tracers 
    9292      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step 
    93       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pwn   ! 3 ocean transport components 
     93      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pu, pv, pw   ! 3 ocean transport components 
    9494      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt_lev1, pt_lev2        ! before and now tracer fields 
    9595      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs             ! tracer trend  
     
    200200               DO jj = 1, jpj 
    201201                  DO ji = 1, jpi 
    202                      zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 
    203                      zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 
     202                     zfp_wk = pw(ji,jj,jk) + ABS( pw(ji,jj,jk) ) 
     203                     zfm_wk = pw(ji,jj,jk) - ABS( pw(ji,jj,jk) ) 
    204204                     ztw(ji,jj,jk) = 0.5_wp * (  zfp_wk * pt_lev1(ji,jj,jk,jn) + zfm_wk * pt_lev1(ji,jj,jk-1,jn)  ) * wmask(ji,jj,jk) 
    205205                  END DO 
     
    210210                  DO jj = 1, jpj 
    211211                     DO ji = 1, jpi 
    212                         ztw(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * pt_lev1(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
     212                        ztw(ji,jj, mikt(ji,jj) ) = pw(ji,jj,mikt(ji,jj)) * pt_lev1(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
    213213                     END DO 
    214214                  END DO    
    215215               ELSE                                ! no cavities: only at the ocean surface 
    216                   ztw(:,:,1) = pwn(:,:,1) * pt_lev1(:,:,1,jn) 
     216                  ztw(:,:,1) = pw(:,:,1) * pt_lev1(:,:,1,jn) 
    217217               ENDIF 
    218218            ENDIF 
     
    233233               DO jj = 1, jpj 
    234234                  DO ji = 1, jpi 
    235                      ztw(ji,jj,jk) = (   0.5_wp * pwn(ji,jj,jk) * ( pt_lev2(ji,jj,jk,jn) + pt_lev2(ji,jj,jk-1,jn) )   & 
     235                     ztw(ji,jj,jk) = (   0.5_wp * pw(ji,jj,jk) * ( pt_lev2(ji,jj,jk,jn) + pt_lev2(ji,jj,jk-1,jn) )   & 
    236236                        &              - ztw(ji,jj,jk)   ) * wmask(ji,jj,jk) 
    237237                  END DO 
     
    248248               DO jj = 2, jpjm1 
    249249                  DO ji = fs_2, fs_jpim1 
    250                      ztw(ji,jj,jk) = pwn(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk) 
    251                   END DO 
    252                END DO 
    253             END DO 
    254             IF( ln_linssh )   ztw(:,:, 1 ) = pwn(:,:,1) * pt_lev2(:,:,1,jn)     !!gm ISF & 4th COMPACT doesn't work 
     250                     ztw(ji,jj,jk) = pw(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk) 
     251                  END DO 
     252               END DO 
     253            END DO 
     254            IF( ln_linssh )   ztw(:,:, 1 ) = pw(:,:,1) * pt_lev2(:,:,1,jn)     !!gm ISF & 4th COMPACT doesn't work 
    255255            ! 
    256256         END SELECT 
     
    269269                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    270270                     zltv(ji,jj,jk) = pt_rhs(ji,jj,jk,jn) - zltv(ji,jj,jk)                          & 
    271                         &           + pt_lev2(ji,jj,jk,jn) * (  pwn(ji,jj,jk) - pwn(ji,jj,jk+1)  )   & 
     271                        &           + pt_lev2(ji,jj,jk,jn) * (  pw(ji,jj,jk) - pw(ji,jj,jk+1)  )   & 
    272272                        &                              * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,ktlev) 
    273273                  END DO 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trabbc.F90

    r10425 r10806  
    5151CONTAINS 
    5252 
    53    SUBROUTINE tra_bbc( kt ) 
     53   SUBROUTINE tra_bbc( kt, ktlev, pts_rhs ) 
    5454      !!---------------------------------------------------------------------- 
    5555      !!                  ***  ROUTINE tra_bbc  *** 
     
    7474      !!---------------------------------------------------------------------- 
    7575      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     76      INTEGER, INTENT(in) ::   ktlev  ! time level index for source terms 
     77      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends 
    7678      ! 
    7779      INTEGER  ::   ji, jj    ! dummy loop indices 
     
    8385      IF( l_trdtra )   THEN         ! Save the input temperature trend 
    8486         ALLOCATE( ztrdt(jpi,jpj,jpk) ) 
    85          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
     87         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) 
    8688      ENDIF 
    8789      !                             !  Add the geothermal trend on temperature 
    8890      DO jj = 2, jpjm1 
    8991         DO ji = 2, jpim1 
    90             tsa(ji,jj,mbkt(ji,jj),jp_tem) = tsa(ji,jj,mbkt(ji,jj),jp_tem) + qgh_trd0(ji,jj) / e3t_n(ji,jj,mbkt(ji,jj)) 
     92            pts_rhs(ji,jj,mbkt(ji,jj),jp_tem) = pts_rhs(ji,jj,mbkt(ji,jj),jp_tem) + qgh_trd0(ji,jj) / e3t(ji,jj,mbkt(ji,jj),ktlev) 
    9193         END DO 
    9294      END DO 
    9395      ! 
    94       CALL lbc_lnk( 'trabbc', tsa(:,:,:,jp_tem) , 'T', 1. ) 
     96      CALL lbc_lnk( 'trabbc', pts_rhs(:,:,:,jp_tem) , 'T', 1. ) 
    9597      ! 
    9698      IF( l_trdtra ) THEN        ! Send the trend for diagnostics 
    97          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
     99         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) - ztrdt(:,:,:) 
    98100         CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbc, ztrdt ) 
    99101         DEALLOCATE( ztrdt ) 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trabbl.F90

    r10425 r10806  
    8989 
    9090 
    91    SUBROUTINE tra_bbl( kt ) 
     91   SUBROUTINE tra_bbl( kt, ktlev1, ktlev2, kt2lev, pts_rhs ) 
    9292      !!---------------------------------------------------------------------- 
    9393      !!                  ***  ROUTINE bbl  *** 
     
    101101      !!              is added to the general tracer trend 
    102102      !!---------------------------------------------------------------------- 
    103       INTEGER, INTENT( in ) ::   kt   ! ocean time-step 
     103      INTEGER, INTENT( in ) ::   kt              ! ocean time-step 
     104      INTEGER, INTENT( in ) ::   ktlev1, ktlev2  ! time level indices for 3-time-level source terms 
     105      INTEGER, INTENT( in ) ::   kt2lev          ! time level index for 2-time-level source terms 
     106      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends 
    104107      ! 
    105108      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds 
     
    110113      IF( l_trdtra )   THEN                         !* Save the T-S input trends 
    111114         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 
    112          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
    113          ztrds(:,:,:) = tsa(:,:,:,jp_sal) 
    114       ENDIF 
    115  
    116       IF( l_bbl )   CALL bbl( kt, nit000, 'TRA' )   !* bbl coef. and transport (only if not already done in trcbbl) 
     115         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) 
     116         ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal) 
     117      ENDIF 
     118 
     119      IF( l_bbl )   CALL bbl( kt, nit000, ktlev1, ktlev2, kt2lev, 'TRA' )   !* bbl coef. and transport (only if not already done in trcbbl) 
    117120 
    118121      IF( nn_bbl_ldf == 1 ) THEN                    !* Diffusive bbl 
    119122         ! 
    120          CALL tra_bbl_dif( tsb, tsa, jpts ) 
     123         CALL tra_bbl_dif( ts(:,:,:,:,ktlev1), e3t(:,:,:,ktlev2), pts_rhs, jpts ) 
    121124         IF( ln_ctl )  & 
    122125         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_ldf  - Ta: ', mask1=tmask, & 
     
    131134      IF( nn_bbl_adv /= 0 ) THEN                    !* Advective bbl 
    132135         ! 
    133          CALL tra_bbl_adv( tsb, tsa, jpts ) 
     136         CALL tra_bbl_adv( ts(:,:,:,:,ktlev1), e3t(:,:,:,ktlev2), pts_rhs, jpts ) 
    134137         IF(ln_ctl)   & 
    135138         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_adv  - Ta: ', mask1=tmask,   & 
     
    143146 
    144147      IF( l_trdtra )   THEN                      ! send the trends for further diagnostics 
    145          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    146          ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
     148         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) - ztrdt(:,:,:) 
     149         ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal) - ztrds(:,:,:) 
    147150         CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbl, ztrdt ) 
    148151         CALL trd_tra( kt, 'TRA', jp_sal, jptra_bbl, ztrds ) 
     
    155158 
    156159 
    157    SUBROUTINE tra_bbl_dif( ptb, pta, kjpt ) 
     160   SUBROUTINE tra_bbl_dif( pt, pe3t, pt_rhs, kjpt ) 
    158161      !!---------------------------------------------------------------------- 
    159162      !!                  ***  ROUTINE tra_bbl_dif  *** 
     
    171174      !!      convection is satified) 
    172175      !! 
    173       !! ** Action  :   pta   increased by the bbl diffusive trend 
     176      !! ** Action  :   pt_rhs   increased by the bbl diffusive trend 
    174177      !! 
    175178      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591. 
     
    177180      !!---------------------------------------------------------------------- 
    178181      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers 
    179       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields 
    180       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta    ! tracer trend 
     182      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt     ! tracer fields 
     183      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pe3t   ! thickness fields 
     184      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs ! tracer trend 
    181185      ! 
    182186      INTEGER  ::   ji, jj, jn   ! dummy loop indices 
     
    191195            DO ji = 1, jpi 
    192196               ik = mbkt(ji,jj)                             ! bottom T-level index 
    193                zptb(ji,jj) = ptb(ji,jj,ik,jn)               ! bottom before T and S 
     197               zptb(ji,jj) = pt(ji,jj,ik,jn)               ! bottom before T and S 
    194198            END DO 
    195199         END DO 
     
    198202            DO ji = 2, jpim1 
    199203               ik = mbkt(ji,jj)                            ! bottom T-level index 
    200                pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn)                                                  & 
     204               pt_rhs(ji,jj,ik,jn) = pt_rhs(ji,jj,ik,jn)                                                  & 
    201205                  &             + (  ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )     & 
    202206                  &                - ahu_bbl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )     & 
    203207                  &                + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )     & 
    204208                  &                - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )  )  & 
    205                   &             * r1_e1e2t(ji,jj) / e3t_n(ji,jj,ik) 
     209                  &             * r1_e1e2t(ji,jj) / pe3t(ji,jj,ik) 
    206210            END DO 
    207211         END DO 
     
    212216 
    213217 
    214    SUBROUTINE tra_bbl_adv( ptb, pta, kjpt ) 
     218   SUBROUTINE tra_bbl_adv( pt, pe3t, pt_rhs, kjpt ) 
    215219      !!---------------------------------------------------------------------- 
    216220      !!                  ***  ROUTINE trc_bbl  *** 
     
    228232      !!---------------------------------------------------------------------- 
    229233      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers 
    230       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields 
    231       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta    ! tracer trend 
     234      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt    ! before and now tracer fields 
     235      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pe3t   ! thickness fields 
     236      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs    ! tracer trend 
    232237      ! 
    233238      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices 
     
    250255                  ! 
    251256                  !                                               ! up  -slope T-point (shelf bottom point) 
    252                   zbtr = r1_e1e2t(iis,jj) / e3t_n(iis,jj,ikus) 
    253                   ztra = zu_bbl * ( ptb(iid,jj,ikus,jn) - ptb(iis,jj,ikus,jn) ) * zbtr 
    254                   pta(iis,jj,ikus,jn) = pta(iis,jj,ikus,jn) + ztra 
     257                  zbtr = r1_e1e2t(iis,jj) / pe3t(iis,jj,ikus) 
     258                  ztra = zu_bbl * ( pt(iid,jj,ikus,jn) - pt(iis,jj,ikus,jn) ) * zbtr 
     259                  pt_rhs(iis,jj,ikus,jn) = pt_rhs(iis,jj,ikus,jn) + ztra 
    255260                  ! 
    256261                  DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column) 
    257                      zbtr = r1_e1e2t(iid,jj) / e3t_n(iid,jj,jk) 
    258                      ztra = zu_bbl * ( ptb(iid,jj,jk+1,jn) - ptb(iid,jj,jk,jn) ) * zbtr 
    259                      pta(iid,jj,jk,jn) = pta(iid,jj,jk,jn) + ztra 
     262                     zbtr = r1_e1e2t(iid,jj) / pe3t(iid,jj,jk) 
     263                     ztra = zu_bbl * ( pt(iid,jj,jk+1,jn) - pt(iid,jj,jk,jn) ) * zbtr 
     264                     pt_rhs(iid,jj,jk,jn) = pt_rhs(iid,jj,jk,jn) + ztra 
    260265                  END DO 
    261266                  ! 
    262                   zbtr = r1_e1e2t(iid,jj) / e3t_n(iid,jj,ikud) 
    263                   ztra = zu_bbl * ( ptb(iis,jj,ikus,jn) - ptb(iid,jj,ikud,jn) ) * zbtr 
    264                   pta(iid,jj,ikud,jn) = pta(iid,jj,ikud,jn) + ztra 
     267                  zbtr = r1_e1e2t(iid,jj) / pe3t(iid,jj,ikud) 
     268                  ztra = zu_bbl * ( pt(iis,jj,ikus,jn) - pt(iid,jj,ikud,jn) ) * zbtr 
     269                  pt_rhs(iid,jj,ikud,jn) = pt_rhs(iid,jj,ikud,jn) + ztra 
    265270               ENDIF 
    266271               ! 
     
    272277                  ! 
    273278                  ! up  -slope T-point (shelf bottom point) 
    274                   zbtr = r1_e1e2t(ji,ijs) / e3t_n(ji,ijs,ikvs) 
    275                   ztra = zv_bbl * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr 
    276                   pta(ji,ijs,ikvs,jn) = pta(ji,ijs,ikvs,jn) + ztra 
     279                  zbtr = r1_e1e2t(ji,ijs) / pe3t(ji,ijs,ikvs) 
     280                  ztra = zv_bbl * ( pt(ji,ijd,ikvs,jn) - pt(ji,ijs,ikvs,jn) ) * zbtr 
     281                  pt_rhs(ji,ijs,ikvs,jn) = pt_rhs(ji,ijs,ikvs,jn) + ztra 
    277282                  ! 
    278283                  DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column) 
    279                      zbtr = r1_e1e2t(ji,ijd) / e3t_n(ji,ijd,jk) 
    280                      ztra = zv_bbl * ( ptb(ji,ijd,jk+1,jn) - ptb(ji,ijd,jk,jn) ) * zbtr 
    281                      pta(ji,ijd,jk,jn) = pta(ji,ijd,jk,jn)  + ztra 
     284                     zbtr = r1_e1e2t(ji,ijd) / pe3t(ji,ijd,jk) 
     285                     ztra = zv_bbl * ( pt(ji,ijd,jk+1,jn) - pt(ji,ijd,jk,jn) ) * zbtr 
     286                     pt_rhs(ji,ijd,jk,jn) = pt_rhs(ji,ijd,jk,jn)  + ztra 
    282287                  END DO 
    283288                  !                                               ! down-slope T-point (deep bottom point) 
    284                   zbtr = r1_e1e2t(ji,ijd) / e3t_n(ji,ijd,ikvd) 
    285                   ztra = zv_bbl * ( ptb(ji,ijs,ikvs,jn) - ptb(ji,ijd,ikvd,jn) ) * zbtr 
    286                   pta(ji,ijd,ikvd,jn) = pta(ji,ijd,ikvd,jn) + ztra 
     289                  zbtr = r1_e1e2t(ji,ijd) / pe3t(ji,ijd,ikvd) 
     290                  ztra = zv_bbl * ( pt(ji,ijs,ikvs,jn) - pt(ji,ijd,ikvd,jn) ) * zbtr 
     291                  pt_rhs(ji,ijd,ikvd,jn) = pt_rhs(ji,ijd,ikvd,jn) + ztra 
    287292               ENDIF 
    288293            END DO 
     
    295300 
    296301 
    297    SUBROUTINE bbl( kt, kit000, cdtype ) 
     302   SUBROUTINE bbl( kt, kit000, ktlev1, ktlev2, kt2lev, cdtype ) 
    298303      !!---------------------------------------------------------------------- 
    299304      !!                  ***  ROUTINE bbl  *** 
     
    323328      INTEGER         , INTENT(in   ) ::   kt       ! ocean time-step index 
    324329      INTEGER         , INTENT(in   ) ::   kit000   ! first time step index 
     330      INTEGER         , INTENT(in   ) ::   ktlev1, ktlev2  ! time level indices for 3-time-levelsource terms 
     331      INTEGER         , INTENT(in   ) ::   kt2lev          ! time level index for 2-time-level source terms 
    325332      CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator) 
    326333      ! 
     
    344351         DO ji = 1, jpi 
    345352            ik = mbkt(ji,jj)                             ! bottom T-level index 
    346             zts (ji,jj,jp_tem) = tsb(ji,jj,ik,jp_tem)    ! bottom before T and S 
    347             zts (ji,jj,jp_sal) = tsb(ji,jj,ik,jp_sal) 
     353            zts (ji,jj,jp_tem) = ts(ji,jj,ik,jp_tem,ktlev1)    ! bottom before T and S 
     354            zts (ji,jj,jp_sal) = ts(ji,jj,ik,jp_sal,ktlev1) 
    348355            ! 
    349             zdep(ji,jj) = gdept_n(ji,jj,ik)              ! bottom T-level reference depth 
    350             zub (ji,jj) = un(ji,jj,mbku(ji,jj))          ! bottom velocity 
    351             zvb (ji,jj) = vn(ji,jj,mbkv(ji,jj)) 
     356            zdep(ji,jj) = gdept(ji,jj,ik,kt2lev)              ! bottom T-level reference depth 
     357            zub (ji,jj) = uu(ji,jj,mbku(ji,jj),ktlev2)          ! bottom velocity 
     358            zvb (ji,jj) = vv(ji,jj,mbkv(ji,jj),ktlev2) 
    352359         END DO 
    353360      END DO 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tradmp.F90

    r10425 r10806  
    7272 
    7373 
    74    SUBROUTINE tra_dmp( kt ) 
     74   SUBROUTINE tra_dmp( kt, ktlev, kt2lev, pts_rhs ) 
    7575      !!---------------------------------------------------------------------- 
    7676      !!                   ***  ROUTINE tra_dmp  *** 
     
    9090      !! ** Action  : - tsa: tracer trends updated with the damping trend 
    9191      !!---------------------------------------------------------------------- 
    92       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     92      INTEGER, INTENT(in) ::   kt     ! ocean time-step index 
     93      INTEGER, INTENT(in) ::   ktlev  ! time level index for 3-time-level source terms 
     94      INTEGER, INTENT(in) ::   kt2lev ! time level index for 2-time-level source terms 
     95      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends 
    9396      ! 
    9497      INTEGER ::   ji, jj, jk, jn   ! dummy loop indices 
     
    101104      IF( l_trdtra )   THEN                    !* Save ta and sa trends 
    102105         ALLOCATE( ztrdts(jpi,jpj,jpk,jpts) )  
    103          ztrdts(:,:,:,:) = tsa(:,:,:,:)  
     106         ztrdts(:,:,:,:) = pts_rhs(:,:,:,:)  
    104107      ENDIF 
    105108      !                           !==  input T-S data at kt  ==! 
     
    113116               DO jj = 2, jpjm1 
    114117                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    115                      tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jn) - tsb(ji,jj,jk,jn) ) 
     118                     pts_rhs(ji,jj,jk,jn) = pts_rhs(ji,jj,jk,jn) + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jn) - ts(ji,jj,jk,jn,ktlev) ) 
    116119                  END DO 
    117120               END DO 
     
    124127               DO ji = fs_2, fs_jpim1   ! vector opt. 
    125128                  IF( avt(ji,jj,jk) <= avt_c ) THEN 
    126                      tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   & 
    127                         &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) 
    128                      tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   & 
    129                         &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 
     129                     pts_rhs(ji,jj,jk,jp_tem) = pts_rhs(ji,jj,jk,jp_tem)   & 
     130                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - ts(ji,jj,jk,jp_tem,ktlev) ) 
     131                     pts_rhs(ji,jj,jk,jp_sal) = pts_rhs(ji,jj,jk,jp_sal)   & 
     132                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - ts(ji,jj,jk,jp_sal,ktlev) ) 
    130133                  ENDIF 
    131134               END DO 
     
    137140            DO jj = 2, jpjm1 
    138141               DO ji = fs_2, fs_jpim1   ! vector opt. 
    139                   IF( gdept_n(ji,jj,jk) >= hmlp (ji,jj) ) THEN 
    140                      tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   & 
    141                         &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) 
    142                      tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   & 
    143                         &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 
     142                  IF( gdept(ji,jj,jk,kt2lev) >= hmlp (ji,jj) ) THEN 
     143                     pts_rhs(ji,jj,jk,jp_tem) = pts_rhs(ji,jj,jk,jp_tem)   & 
     144                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - ts(ji,jj,jk,jp_tem,ktlev) ) 
     145                     pts_rhs(ji,jj,jk,jp_sal) = pts_rhs(ji,jj,jk,jp_sal)   & 
     146                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - ts(ji,jj,jk,jp_sal,ktlev) ) 
    144147                  ENDIF 
    145148               END DO 
     
    150153      ! 
    151154      IF( l_trdtra )   THEN       ! trend diagnostic 
    152          ztrdts(:,:,:,:) = tsa(:,:,:,:) - ztrdts(:,:,:,:) 
     155         ztrdts(:,:,:,:) = pts_rhs(:,:,:,:) - ztrdts(:,:,:,:) 
    153156         CALL trd_tra( kt, 'TRA', jp_tem, jptra_dmp, ztrdts(:,:,:,jp_tem) ) 
    154157         CALL trd_tra( kt, 'TRA', jp_sal, jptra_dmp, ztrdts(:,:,:,jp_sal) ) 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf.F90

    r10068 r10806  
    4747CONTAINS 
    4848 
    49    SUBROUTINE tra_ldf( kt ) 
     49   SUBROUTINE tra_ldf( kt, ktlev1, ktlev2, kt2lev, pts_rhs ) 
    5050      !!---------------------------------------------------------------------- 
    5151      !!                  ***  ROUTINE tra_ldf  *** 
     
    5353      !! ** Purpose :   compute the lateral ocean tracer physics. 
    5454      !!---------------------------------------------------------------------- 
    55       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
     55      INTEGER, INTENT( in ) ::   kt              ! ocean time-step index 
     56      INTEGER, INTENT( in ) ::   ktlev1, ktlev2  ! time level indices for 3-time-level source terms 
     57      INTEGER, INTENT( in ) ::   kt2lev          ! time level index for 2-time-level source terms 
     58      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends 
    5659      !! 
    5760      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds 
     
    6265      IF( l_trdtra )   THEN                    !* Save ta and sa trends 
    6366         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )  
    64          ztrdt(:,:,:) = tsa(:,:,:,jp_tem)  
    65          ztrds(:,:,:) = tsa(:,:,:,jp_sal) 
     67         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem)  
     68         ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal) 
    6669      ENDIF 
    6770      ! 
    6871      SELECT CASE ( nldf_tra )                 !* compute lateral mixing trend and add it to the general trend 
    6972      CASE ( np_lap   )                                  ! laplacian: iso-level operator 
    70          CALL tra_ldf_lap  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb,      tsa, jpts,  1   ) 
     73         CALL tra_ldf_lap  ( kt, nit000, ktlev2, kt2lev, 'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts(:,:,:,:,ktlev1),      pts_rhs, jpts,  1   ) 
    7174      CASE ( np_lap_i )                                  ! laplacian: standard iso-neutral operator (Madec) 
    72          CALL tra_ldf_iso  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsb, tsa, jpts,  1   ) 
     75         CALL tra_ldf_iso  ( kt, nit000, ktlev2, kt2lev, 'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts(:,:,:,:,ktlev1), ts(:,:,:,:,ktlev1), pts_rhs, jpts,  1   ) 
    7376      CASE ( np_lap_it )                                 ! laplacian: triad iso-neutral operator (griffies) 
    74          CALL tra_ldf_triad( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsb, tsa, jpts,  1   ) 
     77         CALL tra_ldf_triad( kt, nit000, ktlev2, kt2lev, 'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts(:,:,:,:,ktlev1), ts(:,:,:,:,ktlev1), pts_rhs, jpts,  1   ) 
    7578      CASE ( np_blp , np_blp_i , np_blp_it )             ! bilaplacian: iso-level & iso-neutral operators 
    76          CALL tra_ldf_blp  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb      , tsa, jpts, nldf_tra ) 
     79         CALL tra_ldf_blp  ( kt, nit000, ktlev2, kt2lev, 'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts(:,:,:,:,ktlev1)      , pts_rhs, jpts, nldf_tra ) 
    7780      END SELECT 
    7881      ! 
    7982      IF( l_trdtra )   THEN                    !* save the horizontal diffusive trends for further diagnostics 
    80          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    81          ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
     83         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) - ztrdt(:,:,:) 
     84         ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal) - ztrds(:,:,:) 
    8285         CALL trd_tra( kt, 'TRA', jp_tem, jptra_ldf, ztrdt ) 
    8386         CALL trd_tra( kt, 'TRA', jp_sal, jptra_ldf, ztrds ) 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_iso.F90

    r10068 r10806  
    4848CONTAINS 
    4949 
    50   SUBROUTINE tra_ldf_iso( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   & 
     50  SUBROUTINE tra_ldf_iso( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, pgu , pgv ,   & 
    5151      &                                                   pgui, pgvi,   & 
    52       &                                       ptb , ptbb, pta , kjpt, kpass ) 
     52      &                                       pt , pt_lev0, pt_rhs , kjpt, kpass ) 
    5353      !!---------------------------------------------------------------------- 
    5454      !!                  ***  ROUTINE tra_ldf_iso  *** 
     
    8787      !!         difft = 1/(e1e2t*e3t) dk[ zftw ] 
    8888      !!      Add this trend to the general trend (ta,sa): 
    89       !!         pta = pta + difft 
    90       !! 
    91       !! ** Action :   Update pta arrays with the before rotated diffusion 
     89      !!         pt_rhs = pt_rhs + difft 
     90      !! 
     91      !! ** Action :   Update pt_rhs arrays with the before rotated diffusion 
    9292      !!---------------------------------------------------------------------- 
    9393      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index 
    9494      INTEGER                              , INTENT(in   ) ::   kit000     ! first time step index 
     95      INTEGER                              , INTENT(in   ) ::   ktlev      ! time level index for e3t 
     96      INTEGER                              , INTENT(in   ) ::   kt2lev     ! time level index for 2-time-level thicknesses 
    9597      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
    9698      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers 
     
    99101      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels 
    100102      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels 
    101       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! tracer (kpass=1) or laplacian of tracer (kpass=2) 
    102       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptbb       ! tracer (only used in kpass=2) 
    103       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend 
     103      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt        ! tracer (kpass=1) or laplacian of tracer (kpass=2) 
     104      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt_lev0       ! tracer (only used in kpass=2) 
     105      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs        ! tracer trend 
    104106      ! 
    105107      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices 
     
    182184                     DO ji = 1, fs_jpim1 
    183185                        akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   & 
    184                            &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) )  ) 
     186                           &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w(ji,jj,jk,kt2lev) * e3w(ji,jj,jk,kt2lev) )  ) 
    185187                     END DO 
    186188                  END DO 
     
    190192                  DO jj = 1, jpjm1 
    191193                     DO ji = 1, fs_jpim1 
    192                         ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) 
     194                        ze3w_2 = e3w(ji,jj,jk,kt2lev) * e3w(ji,jj,jk,kt2lev) 
    193195                        zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  ) 
    194196                        akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt 
     
    219221            DO jj = 1, jpjm1 
    220222               DO ji = 1, fs_jpim1   ! vector opt. 
    221                   zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
    222                   zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 
     223                  zdit(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
     224                  zdjt(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 
    223225               END DO 
    224226            END DO 
     
    248250            ! 
    249251            !                             !== Vertical tracer gradient 
    250             zdk1t(:,:) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * wmask(:,:,jk+1)     ! level jk+1 
     252            zdk1t(:,:) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * wmask(:,:,jk+1)     ! level jk+1 
    251253            ! 
    252254            IF( jk == 1 ) THEN   ;   zdkt(:,:) = zdk1t(:,:)                          ! surface: zdkt(jk=1)=zdkt(jk=2) 
    253             ELSE                 ;   zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk) 
     255            ELSE                 ;   zdkt(:,:) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * wmask(:,:,jk) 
    254256            ENDIF 
    255257            DO jj = 1 , jpjm1            !==  Horizontal fluxes 
    256258               DO ji = 1, fs_jpim1   ! vector opt. 
    257                   zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk) 
    258                   zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk) 
     259                  zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,ktlev) 
     260                  zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,ktlev) 
    259261                  ! 
    260262                  zmsku = 1. / MAX(  wmask(ji+1,jj,jk  ) + wmask(ji,jj,jk+1)   & 
     
    276278            END DO 
    277279            ! 
    278             DO jj = 2 , jpjm1          !== horizontal divergence and add to pta 
     280            DO jj = 2 , jpjm1          !== horizontal divergence and add to pt_rhs 
    279281               DO ji = fs_2, fs_jpim1   ! vector opt. 
    280                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk)      & 
     282                  pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk)      & 
    281283                     &                                           + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )   & 
    282                      &                                        * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     284                     &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,ktlev) 
    283285               END DO 
    284286            END DO 
     
    325327               DO jj = 1, jpjm1 
    326328                  DO ji = fs_2, fs_jpim1 
    327                      ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)   & 
     329                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,kt2lev) * wmask(ji,jj,jk)   & 
    328330                        &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             & 
    329                         &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) 
     331                        &                            * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) 
    330332                  END DO 
    331333               END DO 
     
    340342                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk)    & 
    341343                           &           + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj)   & 
    342                            &           * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) 
     344                           &           * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) / e3w(ji,jj,jk,kt2lev) * wmask(ji,jj,jk) 
    343345                     END DO 
    344346                  END DO 
    345347               END DO  
    346             CASE(  2  )                         ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb  and ptbb gradients, resp. 
     348            CASE(  2  )                         ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt  and pt_lev0 gradients, resp. 
    347349               DO jk = 2, jpkm1  
    348350                  DO jj = 1, jpjm1 
    349351                     DO ji = fs_2, fs_jpim1 
    350                         ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)                      & 
    351                            &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   & 
    352                            &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   ) 
     352                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,kt2lev) * wmask(ji,jj,jk)                      & 
     353                           &                            * (  ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) )   & 
     354                           &                               + akz     (ji,jj,jk) * ( pt_lev0(ji,jj,jk-1,jn) - pt_lev0(ji,jj,jk,jn) )   ) 
    353355                     END DO 
    354356                  END DO 
     
    357359         ENDIF 
    358360         !          
    359          DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pta  ==! 
     361         DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pt_rhs  ==! 
    360362            DO jj = 2, jpjm1 
    361363               DO ji = fs_2, fs_jpim1   ! vector opt. 
    362                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  )   & 
    363                      &                                        * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     364                  pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  )   & 
     365                     &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,ktlev) 
    364366               END DO 
    365367            END DO 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_lap_blp.F90

    r10425 r10806  
    4545CONTAINS 
    4646 
    47    SUBROUTINE tra_ldf_lap( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   & 
     47   SUBROUTINE tra_ldf_lap( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, pgu , pgv ,   & 
    4848      &                                                    pgui, pgvi,   & 
    49       &                                        ptb , pta , kjpt, kpass )  
     49      &                                        pt , pt_rhs , kjpt, kpass )  
    5050      !!---------------------------------------------------------------------- 
    5151      !!                  ***  ROUTINE tra_ldf_lap  *** 
     
    5959      !!          difft = 1/(e1e2t*e3t) {  di-1[ pahu e2u*e3u/e1u di(tb) ] 
    6060      !!                                 + dj-1[ pahv e1v*e3v/e2v dj(tb) ] } 
    61       !!      Add this trend to the general tracer trend pta : 
    62       !!          pta = pta + difft 
    63       !! 
    64       !! ** Action  : - Update pta arrays with the before iso-level  
     61      !!      Add this trend to the general tracer trend pt_rhs : 
     62      !!          pt_rhs = pt_rhs + difft 
     63      !! 
     64      !! ** Action  : - Update pt_rhs arrays with the before iso-level  
    6565      !!                harmonic mixing trend. 
    6666      !!---------------------------------------------------------------------- 
    6767      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index 
    6868      INTEGER                              , INTENT(in   ) ::   kit000     ! first time step index 
     69      INTEGER                              , INTENT(in   ) ::   ktlev      ! time level index for e3t 
     70      INTEGER                              , INTENT(in   ) ::   kt2lev     ! time level index for 2-time-level thicknesses 
    6971      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
    7072      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers 
     
    7375      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels 
    7476      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels 
    75       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! before and now tracer fields 
    76       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend  
     77      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt        ! before and now tracer fields 
     78      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs        ! tracer trend  
    7779      ! 
    7880      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
     
    100102         DO jj = 1, jpjm1 
    101103            DO ji = 1, fs_jpim1   ! vector opt. 
    102                zaheeu(ji,jj,jk) = zsign * pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk)   !!gm   * umask(ji,jj,jk) pah masked! 
    103                zaheev(ji,jj,jk) = zsign * pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk)   !!gm   * vmask(ji,jj,jk) 
     104               zaheeu(ji,jj,jk) = zsign * pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,ktlev)   !!gm   * umask(ji,jj,jk) pah masked! 
     105               zaheev(ji,jj,jk) = zsign * pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,ktlev)   !!gm   * vmask(ji,jj,jk) 
    104106            END DO 
    105107         END DO 
     
    113115            DO jj = 1, jpjm1 
    114116               DO ji = 1, fs_jpim1 
    115                   ztu(ji,jj,jk) = zaheeu(ji,jj,jk) * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) 
    116                   ztv(ji,jj,jk) = zaheev(ji,jj,jk) * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) 
     117                  ztu(ji,jj,jk) = zaheeu(ji,jj,jk) * ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) 
     118                  ztv(ji,jj,jk) = zaheev(ji,jj,jk) * ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) 
    117119               END DO 
    118120            END DO 
     
    138140            DO jj = 2, jpjm1 
    139141               DO ji = fs_2, fs_jpim1 
    140                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)     & 
     142                  pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)     & 
    141143                     &                                   + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) )   & 
    142                      &                                / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) 
     144                     &                                / ( e1e2t(ji,jj) * e3t(ji,jj,jk,ktlev) ) 
    143145               END DO 
    144146            END DO 
     
    159161    
    160162 
    161    SUBROUTINE tra_ldf_blp( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   & 
     163   SUBROUTINE tra_ldf_blp( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, pgu , pgv ,   & 
    162164      &                                                    pgui, pgvi,   & 
    163       &                                                    ptb , pta , kjpt, kldf ) 
     165      &                                                    pt , pt_rhs , kjpt, kldf ) 
    164166      !!---------------------------------------------------------------------- 
    165167      !!                 ***  ROUTINE tra_ldf_blp  *** 
     
    172174      !!      It is computed by two successive calls to laplacian routine 
    173175      !! 
    174       !! ** Action :   pta   updated with the before rotated bilaplacian diffusion 
     176      !! ** Action :   pt_rhs   updated with the before rotated bilaplacian diffusion 
    175177      !!---------------------------------------------------------------------- 
    176178      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index 
    177179      INTEGER                              , INTENT(in   ) ::   kit000     ! first time step index 
     180      INTEGER                              , INTENT(in   ) ::   ktlev      ! time level index for e3t 
     181      INTEGER                              , INTENT(in   ) ::   kt2lev     ! time level index for 2-time-level thicknesses 
    178182      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
    179183      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers 
     
    182186      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels 
    183187      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top levels 
    184       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! before and now tracer fields 
    185       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend 
     188      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt        ! before and now tracer fields 
     189      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs        ! tracer trend 
    186190      ! 
    187191      INTEGER ::   ji, jj, jk, jn   ! dummy loop indices 
     
    203207      zlap(:,:,:,:) = 0._wp 
    204208      ! 
    205       SELECT CASE ( kldf )       !==  1st laplacian applied to ptb (output in zlap)  ==! 
     209      SELECT CASE ( kldf )       !==  1st laplacian applied to pt (output in zlap)  ==! 
    206210      ! 
    207211      CASE ( np_blp    )               ! iso-level bilaplacian 
    208          CALL tra_ldf_lap  ( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb,      zlap, kjpt, 1 ) 
     212         CALL tra_ldf_lap  ( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, pt,      zlap, kjpt, 1 ) 
    209213      CASE ( np_blp_i  )               ! rotated   bilaplacian : standard operator (Madec) 
    210          CALL tra_ldf_iso  ( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, ptb, zlap, kjpt, 1 ) 
     214         CALL tra_ldf_iso  ( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, pt, pt, zlap, kjpt, 1 ) 
    211215      CASE ( np_blp_it )               ! rotated  bilaplacian : triad operator (griffies) 
    212          CALL tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, ptb, zlap, kjpt, 1 ) 
     216         CALL tra_ldf_triad( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, pt, pt, zlap, kjpt, 1 ) 
    213217      END SELECT 
    214218      ! 
     
    219223      ENDIF 
    220224      ! 
    221       SELECT CASE ( kldf )       !==  2nd laplacian applied to zlap (output in pta)  ==! 
     225      SELECT CASE ( kldf )       !==  2nd laplacian applied to zlap (output in pt_rhs)  ==! 
    222226      ! 
    223227      CASE ( np_blp    )               ! iso-level bilaplacian 
    224          CALL tra_ldf_lap  ( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pta,      kjpt, 2 ) 
     228         CALL tra_ldf_lap  ( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pt_rhs,      kjpt, 2 ) 
    225229      CASE ( np_blp_i  )               ! rotated   bilaplacian : standard operator (Madec) 
    226          CALL tra_ldf_iso  ( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, ptb, pta, kjpt, 2 ) 
     230         CALL tra_ldf_iso  ( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pt, pt_rhs, kjpt, 2 ) 
    227231      CASE ( np_blp_it )               ! rotated  bilaplacian : triad operator (griffies) 
    228          CALL tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, ptb, pta, kjpt, 2 ) 
     232         CALL tra_ldf_triad( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pt, pt_rhs, kjpt, 2 ) 
    229233      END SELECT 
    230234      ! 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_triad.F90

    r10425 r10806  
    4848CONTAINS 
    4949 
    50   SUBROUTINE tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   & 
     50  SUBROUTINE tra_ldf_triad( kt, kit000, ktlev, kt2lev, cdtype, pahu, pahv, pgu , pgv ,   & 
    5151      &                                                     pgui, pgvi,   & 
    52       &                                         ptb , ptbb, pta , kjpt, kpass ) 
     52      &                                         pt , pt_lev0, pt_rhs , kjpt, kpass ) 
    5353      !!---------------------------------------------------------------------- 
    5454      !!                  ***  ROUTINE tra_ldf_triad  *** 
     
    6666      !!      see documentation for the desciption 
    6767      !! 
    68       !! ** Action :   pta   updated with the before rotated diffusion 
     68      !! ** Action :   pt_rhs   updated with the before rotated diffusion 
    6969      !!               ah_wslp2 .... 
    7070      !!               akz   stabilizing vertical diffusivity coefficient (used in trazdf_imp) 
     
    7272      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index 
    7373      INTEGER                              , INTENT(in   ) ::   kit000     ! first time step index 
     74      INTEGER                              , INTENT(in   ) ::   ktlev      ! time level index for e3t 
     75      INTEGER                              , INTENT(in   ) ::   kt2lev     ! time level index for 2-time-level thicknesses 
    7476      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
    7577      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers 
     
    7880      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu , pgv  ! tracer gradient at pstep levels 
    7981      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels 
    80       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! tracer (kpass=1) or laplacian of tracer (kpass=2) 
    81       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptbb       ! tracer (only used in kpass=2) 
    82       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend 
     82      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt        ! tracer (kpass=1) or laplacian of tracer (kpass=2) 
     83      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt_lev0       ! tracer (only used in kpass=2) 
     84      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs        ! tracer trend 
    8385      ! 
    8486      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices 
     
    142144                  DO jj = 1, jpjm1 
    143145                     DO ji = 1, fs_jpim1 
    144                         ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp) 
    145                         zbu   = e1e2u(ji,jj) * e3u_n(ji,jj,jk) 
     146                        ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,kt2lev) 
     147                        zbu   = e1e2u(ji,jj) * e3u(ji,jj,jk,ktlev) 
    146148                        zah   = 0.25_wp * pahu(ji,jj,jk) 
    147149                        zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 
    148150                        ! Subtract s-coordinate slope at t-points to give slope rel to s-surfaces (do this by *adding* gradient of depth) 
    149                         zslope2 = zslope_skew + ( gdept_n(ji+1,jj,jk) - gdept_n(ji,jj,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp) 
     151                        zslope2 = zslope_skew + ( gdept(ji+1,jj,jk,kt2lev) - gdept(ji,jj,jk,kt2lev) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp) 
    150152                        zslope2 = zslope2 *zslope2 
    151153                        ah_wslp2(ji+ip,jj,jk+kp) = ah_wslp2(ji+ip,jj,jk+kp) + zah * zbu * ze3wr * r1_e1e2t(ji+ip,jj) * zslope2 
     
    166168                  DO jj = 1, jpjm1 
    167169                     DO ji = 1, fs_jpim1 
    168                         ze3wr = 1.0_wp / e3w_n(ji,jj+jp,jk+kp) 
    169                         zbv   = e1e2v(ji,jj) * e3v_n(ji,jj,jk) 
     170                        ze3wr = 1.0_wp / e3w(ji,jj+jp,jk+kp,kt2lev) 
     171                        zbv   = e1e2v(ji,jj) * e3v(ji,jj,jk,ktlev) 
    170172                        zah   = 0.25_wp * pahv(ji,jj,jk) 
    171173                        zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 
    172174                        ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces 
    173175                        !    (do this by *adding* gradient of depth) 
    174                         zslope2 = zslope_skew + ( gdept_n(ji,jj+1,jk) - gdept_n(ji,jj,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp) 
     176                        zslope2 = zslope_skew + ( gdept(ji,jj+1,jk,kt2lev) - gdept(ji,jj,jk,kt2lev) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp) 
    175177                        zslope2 = zslope2 * zslope2 
    176178                        ah_wslp2(ji,jj+jp,jk+kp) = ah_wslp2(ji,jj+jp,jk+kp) + zah * zbv * ze3wr * r1_e1e2t(ji,jj+jp) * zslope2 
     
    193195                     DO ji = 1, fs_jpim1 
    194196                        akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   & 
    195                            &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) )  ) 
     197                           &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w(ji,jj,jk,kt2lev) * e3w(ji,jj,jk,kt2lev) )  ) 
    196198                     END DO 
    197199                  END DO 
     
    201203                  DO jj = 1, jpjm1 
    202204                     DO ji = 1, fs_jpim1 
    203                         ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) 
     205                        ze3w_2 = e3w(ji,jj,jk,kt2lev) * e3w(ji,jj,jk,kt2lev) 
    204206                        zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  ) 
    205207                        akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt 
     
    229231            DO jj = 1, jpjm1 
    230232               DO ji = 1, fs_jpim1   ! vector opt. 
    231                   zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
    232                   zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 
     233                  zdit(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
     234                  zdjt(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 
    233235               END DO 
    234236            END DO 
     
    257259         DO jk = 1, jpkm1 
    258260            !                    !==  Vertical tracer gradient at level jk and jk+1 
    259             zdkt3d(:,:,1) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * tmask(:,:,jk+1) 
     261            zdkt3d(:,:,1) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * tmask(:,:,jk+1) 
    260262            ! 
    261263            !                    ! surface boundary condition: zdkt3d(jk=0)=zdkt3d(jk=1) 
    262264            IF( jk == 1 ) THEN   ;   zdkt3d(:,:,0) = zdkt3d(:,:,1) 
    263             ELSE                 ;   zdkt3d(:,:,0) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * tmask(:,:,jk) 
     265            ELSE                 ;   zdkt3d(:,:,0) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * tmask(:,:,jk) 
    264266            ENDIF 
    265267            ! 
     
    273275                           ze1ur = r1_e1u(ji,jj) 
    274276                           zdxt  = zdit(ji,jj,jk) * ze1ur 
    275                            ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp) 
     277                           ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,kt2lev) 
    276278                           zdzt  = zdkt3d(ji+ip,jj,kp) * ze3wr 
    277279                           zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 
    278280                           zslope_iso  = triadi  (ji+ip,jj,jk,1-ip,kp) 
    279281                           ! 
    280                            zbu = 0.25_wp * e1e2u(ji,jj) * e3u_n(ji,jj,jk) 
     282                           zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,ktlev) 
    281283                           ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????   ahu is masked.... 
    282284                           zah = pahu(ji,jj,jk) 
     
    296298                           ze2vr = r1_e2v(ji,jj) 
    297299                           zdyt  = zdjt(ji,jj,jk) * ze2vr 
    298                            ze3wr = 1._wp / e3w_n(ji,jj+jp,jk+kp) 
     300                           ze3wr = 1._wp / e3w(ji,jj+jp,jk+kp,kt2lev) 
    299301                           zdzt  = zdkt3d(ji,jj+jp,kp) * ze3wr 
    300302                           zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 
    301303                           zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp) 
    302                            zbv = 0.25_wp * e1e2v(ji,jj) * e3v_n(ji,jj,jk) 
     304                           zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,ktlev) 
    303305                           ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????  ahv is masked... 
    304306                           zah = pahv(ji,jj,jk) 
     
    320322                           ze1ur = r1_e1u(ji,jj) 
    321323                           zdxt  = zdit(ji,jj,jk) * ze1ur 
    322                            ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp) 
     324                           ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,kt2lev) 
    323325                           zdzt  = zdkt3d(ji+ip,jj,kp) * ze3wr 
    324326                           zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 
    325327                           zslope_iso  = triadi(ji+ip,jj,jk,1-ip,kp) 
    326328                           ! 
    327                            zbu = 0.25_wp * e1e2u(ji,jj) * e3u_n(ji,jj,jk) 
     329                           zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,ktlev) 
    328330                           ! ln_botmix_triad is .F. mask zah for bottom half cells 
    329331                           zah = pahu(ji,jj,jk) * umask(ji,jj,jk+kp)         ! pahu(ji+ip,jj,jk)   ===>>  ???? 
     
    343345                           ze2vr = r1_e2v(ji,jj) 
    344346                           zdyt  = zdjt(ji,jj,jk) * ze2vr 
    345                            ze3wr = 1._wp / e3w_n(ji,jj+jp,jk+kp) 
     347                           ze3wr = 1._wp / e3w(ji,jj+jp,jk+kp,kt2lev) 
    346348                           zdzt  = zdkt3d(ji,jj+jp,kp) * ze3wr 
    347349                           zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 
    348350                           zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp) 
    349                            zbv = 0.25_wp * e1e2v(ji,jj) * e3v_n(ji,jj,jk) 
     351                           zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,ktlev) 
    350352                           ! ln_botmix_triad is .F. mask zah for bottom half cells 
    351353                           zah = pahv(ji,jj,jk) * vmask(ji,jj,jk+kp)         ! pahv(ji,jj+jp,jk)  ???? 
     
    362364            DO jj = 2 , jpjm1 
    363365               DO ji = fs_2, fs_jpim1  ! vector opt. 
    364                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  zftu(ji-1,jj,jk) - zftu(ji,jj,jk)       & 
     366                  pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  zftu(ji-1,jj,jk) - zftu(ji,jj,jk)       & 
    365367                     &                                           + zftv(ji,jj-1,jk) - zftv(ji,jj,jk)   )   & 
    366                      &                                        / (  e1e2t(ji,jj) * e3t_n(ji,jj,jk)  ) 
     368                     &                                        / (  e1e2t(ji,jj) * e3t(ji,jj,jk,ktlev)  ) 
    367369               END DO 
    368370            END DO 
     
    375377               DO jj = 1, jpjm1 
    376378                  DO ji = fs_2, fs_jpim1 
    377                      ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk)   & 
     379                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,kt2lev) * tmask(ji,jj,jk)   & 
    378380                        &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             & 
    379                         &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) 
     381                        &                            * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) 
    380382                  END DO 
    381383               END DO 
     
    387389                  DO jj = 1, jpjm1 
    388390                     DO ji = fs_2, fs_jpim1 
    389                         ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk)             & 
    390                            &                            * ah_wslp2(ji,jj,jk) * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) 
     391                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,kt2lev) * tmask(ji,jj,jk)             & 
     392                           &                            * ah_wslp2(ji,jj,jk) * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) 
    391393                     END DO 
    392394                  END DO 
    393395               END DO  
    394             CASE(  2  )                            ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb  and ptbb gradients, resp. 
     396            CASE(  2  )                            ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt  and pt_lev0 gradients, resp. 
    395397               DO jk = 2, jpkm1  
    396398                  DO jj = 1, jpjm1 
    397399                     DO ji = fs_2, fs_jpim1 
    398                         ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk)                      & 
    399                            &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   & 
    400                            &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   ) 
     400                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,kt2lev) * tmask(ji,jj,jk)                      & 
     401                           &                            * (  ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) )   & 
     402                           &                               + akz     (ji,jj,jk) * ( pt_lev0(ji,jj,jk-1,jn) - pt_lev0(ji,jj,jk,jn) )   ) 
    401403                     END DO 
    402404                  END DO 
     
    405407         ENDIF 
    406408         ! 
    407          DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pta  ==! 
     409         DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pt_rhs  ==! 
    408410            DO jj = 2, jpjm1 
    409411               DO ji = fs_2, fs_jpim1  ! vector opt. 
    410                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   & 
    411                      &                                        / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) 
     412                  pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   & 
     413                     &                                        / ( e1e2t(ji,jj) * e3t(ji,jj,jk,ktlev) ) 
    412414               END DO 
    413415            END DO 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traqsr.F90

    r10425 r10806  
    7575CONTAINS 
    7676 
    77    SUBROUTINE tra_qsr( kt ) 
     77   SUBROUTINE tra_qsr( kt, ktlev, kt2lev, pts_rhs ) 
    7878      !!---------------------------------------------------------------------- 
    7979      !!                  ***  ROUTINE tra_qsr  *** 
     
    102102      !!---------------------------------------------------------------------- 
    103103      INTEGER, INTENT(in) ::   kt     ! ocean time-step 
     104      INTEGER, INTENT(in) ::   ktlev  ! time level index for 3-time-level source terms 
     105      INTEGER, INTENT(in) ::   kt2lev ! time level index for 2-time-level source terms 
     106      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends 
    104107      ! 
    105108      INTEGER  ::   ji, jj, jk               ! dummy loop indices 
     
    126129      IF( l_trdtra ) THEN      ! trends diagnostic: save the input temperature trend 
    127130         ALLOCATE( ztrdt(jpi,jpj,jpk) )  
    128          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
     131         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) 
    129132      ENDIF 
    130133      ! 
     
    172175                     zze     = 568.2 * zCtot**(-0.746) 
    173176                     IF( zze > 102. ) zze = 200.0 * zCtot**(-0.293) 
    174                      zpsi    = gdepw_n(ji,jj,jk) / zze 
     177                     zpsi    = gdepw(ji,jj,jk,kt2lev) / zze 
    175178                     ! 
    176179                     zlogc   = LOG( zchl ) 
     
    218221            DO jj = 2, jpjm1 
    219222               DO ji = fs_2, fs_jpim1 
    220                   zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * xsi0r       ) 
    221                   zc1 = ze1(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekb(ji,jj) ) 
    222                   zc2 = ze2(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekg(ji,jj) ) 
    223                   zc3 = ze3(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekr(ji,jj) ) 
     223                  zc0 = ze0(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,ktlev) * xsi0r       ) 
     224                  zc1 = ze1(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,ktlev) * zekb(ji,jj) ) 
     225                  zc2 = ze2(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,ktlev) * zekg(ji,jj) ) 
     226                  zc3 = ze3(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,ktlev) * zekr(ji,jj) ) 
    224227                  ze0(ji,jj,jk) = zc0 
    225228                  ze1(ji,jj,jk) = zc1 
     
    248251            DO jj = 2, jpjm1 
    249252               DO ji = fs_2, fs_jpim1 
    250                   zc0 = zz0 * EXP( -gdepw_n(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -gdepw_n(ji,jj,jk  )*xsi1r ) 
    251                   zc1 = zz0 * EXP( -gdepw_n(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -gdepw_n(ji,jj,jk+1)*xsi1r ) 
     253                  zc0 = zz0 * EXP( -gdepw(ji,jj,jk  ,kt2lev)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk  ,kt2lev)*xsi1r ) 
     254                  zc1 = zz0 * EXP( -gdepw(ji,jj,jk+1,kt2lev)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk+1,kt2lev)*xsi1r ) 
    252255                  qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) )  
    253256               END DO 
     
    261264         DO jj = 2, jpjm1        !-----------------------------! 
    262265            DO ji = fs_2, fs_jpim1   ! vector opt. 
    263                tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   & 
    264                   &                 + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t_n(ji,jj,jk) 
     266               pts_rhs(ji,jj,jk,jp_tem) = pts_rhs(ji,jj,jk,jp_tem)   & 
     267                  &                 + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t(ji,jj,jk,ktlev) 
    265268            END DO 
    266269         END DO 
     
    295298      ! 
    296299      IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics 
    297          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
     300         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) - ztrdt(:,:,:) 
    298301         CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrdt ) 
    299302         DEALLOCATE( ztrdt )  
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trasbc.F90

    r10499 r10806  
    5151CONTAINS 
    5252 
    53    SUBROUTINE tra_sbc ( kt ) 
     53   SUBROUTINE tra_sbc ( kt, ktlev, pts_rhs ) 
    5454      !!---------------------------------------------------------------------- 
    5555      !!                  ***  ROUTINE tra_sbc  *** 
     
    6363      !!      (2) Fwe , tracer carried with the water that is exchanged with air+ice.  
    6464      !!               The input forcing fields (emp, rnf, sfx, isf) contain Fext+Fwe, 
    65       !!             they are simply added to the tracer trend (tsa). 
     65      !!             they are simply added to the tracer trend (pts_rhs). 
    6666      !!               In linear free surface case (ln_linssh=T), the volume of the 
    6767      !!             ocean does not change with the water exchanges at the (air+ice)-sea 
     
    6969      !!             concentration/dilution effect associated with water exchanges. 
    7070      !! 
    71       !! ** Action  : - Update tsa with the surface boundary condition trend  
     71      !! ** Action  : - Update pts_rhs with the surface boundary condition trend  
    7272      !!              - send trends to trdtra module for further diagnostics(l_trdtra=T) 
    7373      !!---------------------------------------------------------------------- 
    74       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     74      INTEGER, INTENT(in) ::   kt     ! ocean time-step index 
     75      INTEGER, INTENT(in) ::   ktlev  ! time level index for source terms 
     76      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends 
    7577      ! 
    7678      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices   
     
    9092      IF( l_trdtra ) THEN                    !* Save ta and sa trends 
    9193         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )  
    92          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
    93          ztrds(:,:,:) = tsa(:,:,:,jp_sal) 
     94         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) 
     95         ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal) 
    9496      ENDIF 
    9597      ! 
     
    131133         DO jj = 2, jpj                         !==>> add concentration/dilution effect due to constant volume cell 
    132134            DO ji = fs_2, fs_jpim1   ! vector opt. 
    133                sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rau0 * emp(ji,jj) * tsn(ji,jj,1,jp_tem) 
    134                sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rau0 * emp(ji,jj) * tsn(ji,jj,1,jp_sal) 
     135               sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rau0 * emp(ji,jj) * ts(ji,jj,1,jp_tem,ktlev) 
     136               sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rau0 * emp(ji,jj) * ts(ji,jj,1,jp_sal,ktlev) 
    135137            END DO 
    136138         END DO                                 !==>> output c./d. term 
    137          IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * tsn(:,:,1,jp_tem) ) 
    138          IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * tsn(:,:,1,jp_sal) ) 
     139         IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * ts(:,:,1,jp_tem,ktlev) ) 
     140         IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * ts(:,:,1,jp_sal,ktlev) ) 
    139141      ENDIF 
    140142      ! 
     
    142144         DO jj = 2, jpj 
    143145            DO ji = fs_2, fs_jpim1   ! vector opt.   
    144                tsa(ji,jj,1,jn) = tsa(ji,jj,1,jn) + zfact * ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) / e3t_n(ji,jj,1) 
     146               pts_rhs(ji,jj,1,jn) = pts_rhs(ji,jj,1,jn) + zfact * ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) / e3t(ji,jj,1,ktlev) 
    145147            END DO 
    146148         END DO 
     
    173175               DO jk = ikt, ikb - 1 
    174176               ! compute trend 
    175                   tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                                & 
     177                  pts_rhs(ji,jj,jk,jp_tem) = pts_rhs(ji,jj,jk,jp_tem)                                                & 
    176178                     &           + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) )             & 
    177179                     &           * r1_hisf_tbl(ji,jj) 
     
    180182               ! level partially include in ice shelf boundary layer  
    181183               ! compute trend 
    182                tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                                 & 
     184               pts_rhs(ji,jj,ikb,jp_tem) = pts_rhs(ji,jj,ikb,jp_tem)                                                 & 
    183185                  &              + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) )             & 
    184186                  &              * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) 
     
    199201                  zdep = zfact / h_rnf(ji,jj) 
    200202                  DO jk = 1, nk_rnf(ji,jj) 
    201                                         tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                 & 
     203                                        pts_rhs(ji,jj,jk,jp_tem) = pts_rhs(ji,jj,jk,jp_tem)                                 & 
    202204                                           &                 +  ( rnf_tsc_b(ji,jj,jp_tem) + rnf_tsc(ji,jj,jp_tem) ) * zdep 
    203                      IF( ln_rnf_sal )   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)                                 & 
     205                     IF( ln_rnf_sal )   pts_rhs(ji,jj,jk,jp_sal) = pts_rhs(ji,jj,jk,jp_sal)                                 & 
    204206                                           &                 +  ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep  
    205207                  END DO 
     
    209211      ENDIF 
    210212 
    211       IF( iom_use('rnf_x_sst') )   CALL iom_put( "rnf_x_sst", rnf*tsn(:,:,1,jp_tem) )   ! runoff term on sst 
    212       IF( iom_use('rnf_x_sss') )   CALL iom_put( "rnf_x_sss", rnf*tsn(:,:,1,jp_sal) )   ! runoff term on sss 
     213      IF( iom_use('rnf_x_sst') )   CALL iom_put( "rnf_x_sst", rnf*ts(:,:,1,jp_tem,ktlev) )   ! runoff term on sst 
     214      IF( iom_use('rnf_x_sss') )   CALL iom_put( "rnf_x_sss", rnf*ts(:,:,1,jp_sal,ktlev) )   ! runoff term on sss 
    213215 
    214216#if defined key_asminc 
     
    223225            DO jj = 2, jpj  
    224226               DO ji = fs_2, fs_jpim1 
    225                   ztim = ssh_iau(ji,jj) / e3t_n(ji,jj,1) 
    226                   tsa(ji,jj,1,jp_tem) = tsa(ji,jj,1,jp_tem) + tsn(ji,jj,1,jp_tem) * ztim 
    227                   tsa(ji,jj,1,jp_sal) = tsa(ji,jj,1,jp_sal) + tsn(ji,jj,1,jp_sal) * ztim 
     227                  ztim = ssh_iau(ji,jj) / e3t(ji,jj,1,ktlev) 
     228                  pts_rhs(ji,jj,1,jp_tem) = pts_rhs(ji,jj,1,jp_tem) + ts(ji,jj,1,jp_tem,ktlev) * ztim 
     229                  pts_rhs(ji,jj,1,jp_sal) = pts_rhs(ji,jj,1,jp_sal) + ts(ji,jj,1,jp_sal,ktlev) * ztim 
    228230               END DO 
    229231            END DO 
     
    232234               DO ji = fs_2, fs_jpim1 
    233235                  ztim = ssh_iau(ji,jj) / ( ht_n(ji,jj) + 1. - ssmask(ji, jj) ) 
    234                   tsa(ji,jj,:,jp_tem) = tsa(ji,jj,:,jp_tem) + tsn(ji,jj,:,jp_tem) * ztim 
    235                   tsa(ji,jj,:,jp_sal) = tsa(ji,jj,:,jp_sal) + tsn(ji,jj,:,jp_sal) * ztim 
     236                  pts_rhs(ji,jj,:,jp_tem) = pts_rhs(ji,jj,:,jp_tem) + ts(ji,jj,:,jp_tem,ktlev) * ztim 
     237                  pts_rhs(ji,jj,:,jp_sal) = pts_rhs(ji,jj,:,jp_sal) + ts(ji,jj,:,jp_sal,ktlev) * ztim 
    236238               END DO   
    237239            END DO   
     
    250252            DO jj = 2, jpj  
    251253               DO ji = fs_2, fs_jpim1 
    252                   zdep = 1._wp / e3t_n(ji,jj,jk)  
    253                   tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) - htsc_iscpl(ji,jj,jk,jp_tem) * zdep 
    254                   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) - htsc_iscpl(ji,jj,jk,jp_sal) * zdep   
     254                  zdep = 1._wp / e3t(ji,jj,jk,ktlev)  
     255                  pts_rhs(ji,jj,jk,jp_tem) = pts_rhs(ji,jj,jk,jp_tem) - htsc_iscpl(ji,jj,jk,jp_tem) * zdep 
     256                  pts_rhs(ji,jj,jk,jp_sal) = pts_rhs(ji,jj,jk,jp_sal) - htsc_iscpl(ji,jj,jk,jp_sal) * zdep   
    255257               END DO   
    256258            END DO   
     
    259261 
    260262      IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics 
    261          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    262          ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
     263         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) - ztrdt(:,:,:) 
     264         ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal) - ztrds(:,:,:) 
    263265         CALL trd_tra( kt, 'TRA', jp_tem, jptra_nsr, ztrdt ) 
    264266         CALL trd_tra( kt, 'TRA', jp_sal, jptra_nsr, ztrds ) 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/oce.F90

    r10802 r10806  
    2222   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:),   TARGET ::   ww             !: vertical velocity            [m/s] 
    2323   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::           wi             !: vertical vel. (adaptive-implicit) [m/s] 
    24    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::           hdivn          !: horizontal divergence        [s-1] 
     24   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET   ::   hdiv           !: horizontal divergence        [s-1] 
    2525   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:), TARGET :: ts             !: 4D T-S fields                  [Celsius,psu]  
    26    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   rab_b,  rab_n          !: thermal/haline expansion coef. [Celsius-1,psu-1] 
    27    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   rn2b ,  rn2            !: brunt-vaisala frequency**2     [s-2] 
     26   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:), TARGET :: r_ab           !: thermal/haline expansion coef. [Celsius-1,psu-1] 
     27   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:),TARGET    :: r_n2           !: brunt-vaisala frequency**2     [s-2] 
    2828   ! 
    2929   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rhd    !: in situ density anomalie rhd=(rho-rau0)/rau0  [no units] 
     
    7272   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:)   ::   vb   ,  vn    , va     !: j-horizontal velocity        [m/s] 
    7373   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:)   ::           wn             !: k-vertical   velocity        [m/s] 
     74   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:)   ::           hdivn          !: horizontal divergence        [s-1] 
     75   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:,:) ::   rab_b,  rab_n          !: thermal/haline expansion coef. [Celsius-1,psu-1] 
     76   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:)   ::   rn2b ,  rn2            !: brunt-vaisala frequency**2     [s-2] 
    7477   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:,:) ::   tsb  ,  tsn   , tsa    !: 4D T-S fields                [Celsius,psu]          
    7578   !! TEMPORARY POINTERS FOR DEVELOPMENT ONLY 
     
    9194      ierr(:) = 0  
    9295      ALLOCATE( uu   (jpi,jpj,jpk, jpt) , vv   (jpi,jpj,jpk,jpt)  ,                             & 
    93          &      ww   (jpi,jpj,jpk)      , hdivn(jpi,jpj,jpk)      ,                             & 
     96         &      ww   (jpi,jpj,jpk)      , hdiv(jpi,jpj,jpk)       ,                             & 
    9497         &      ts  (jpi,jpj,jpk,jpts,jpt) ,                                                    & 
    95          &      rab_b(jpi,jpj,jpk,jpts) , rab_n(jpi,jpj,jpk,jpts) ,                             & 
    96          &      rn2b (jpi,jpj,jpk)      , rn2  (jpi,jpj,jpk)      ,                             & 
    97          &      rhd  (jpi,jpj,jpk)      , rhop (jpi,jpj,jpk)                              , STAT=ierr(1) ) 
     98         &      r_ab (jpi,jpj,jpk,jpts,jpt), r_n2 (jpi,jpj,jpk,jpt),                            & 
     99         &      rhd  (jpi,jpj,jpk)         , rhop (jpi,jpj,jpk)                              , STAT=ierr(1) ) 
    98100         ! 
    99101      ALLOCATE( sshb(jpi,jpj)     , sshn(jpi,jpj)   , ssha(jpi,jpj)   ,     & 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90

    r10802 r10806  
    178178!!jc: fs simplification 
    179179                             
    180                          ua(:,:,:) = 0._wp            ! set dynamics trends to zero 
    181                          va(:,:,:) = 0._wp 
     180                         uu(:,:,:,Nrhs) = 0._wp            ! set dynamics trends to zero 
     181                         vv(:,:,:,Nrhs) = 0._wp 
    182182 
    183183      IF(  lk_asminc .AND. ln_asmiau .AND. ln_dyninc )   & 
     
    190190                         CALL dyn_adv       ( kstp, Nm1, Nnn, uu(:,:,:,Nrhs), vv(:,:,:,Nrhs) )  ! advection (vector or flux form) 
    191191                         CALL dyn_vor       ( kstp,      Nnn, uu(:,:,:,Nrhs), vv(:,:,:,Nrhs) )  ! vorticity term including Coriolis 
    192                          CALL dyn_ldf       ( kstp )  ! lateral mixing 
     192                         CALL dyn_ldf       ( kstp, Nm1, Nnn, uu(:,:,:,Nrhs), vv(:,:,:,Nrhs) )  ! lateral mixing 
    193193      IF( ln_zdfosm  )   CALL dyn_osm       ( kstp )  ! OSMOSIS non-local velocity fluxes 
    194194                         CALL dyn_hpg       ( kstp )  ! horizontal gradient of Hydrostatic pressure 
     
    233233      ! Active tracers                               
    234234      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    235                          tsa(:,:,:,:) = 0._wp         ! set tracer trends to zero 
     235                         ts(:,:,:,:,Nrhs) = 0._wp         ! set tracer trends to zero 
    236236 
    237237      IF(  lk_asminc .AND. ln_asmiau .AND. & 
    238238         & ln_trainc )   CALL tra_asm_inc   ( kstp )  ! apply tracer assimilation increment 
    239                          CALL tra_sbc       ( kstp )  ! surface boundary condition 
    240       IF( ln_traqsr  )   CALL tra_qsr       ( kstp )  ! penetrative solar radiation qsr 
    241       IF( ln_trabbc  )   CALL tra_bbc       ( kstp )  ! bottom heat flux 
    242       IF( ln_trabbl  )   CALL tra_bbl       ( kstp )  ! advective (and/or diffusive) bottom boundary layer scheme 
    243       IF( ln_tradmp  )   CALL tra_dmp       ( kstp )  ! internal damping trends 
     239                         CALL tra_sbc       ( kstp, Nnn, ts(:,:,:,:,Nrhs) )  ! surface boundary condition 
     240      IF( ln_traqsr  )   CALL tra_qsr       ( kstp, Nnn, Nnn_2lev, ts(:,:,:,:,Nrhs) )  ! penetrative solar radiation qsr 
     241      IF( ln_trabbc  )   CALL tra_bbc       ( kstp, Nnn, ts(:,:,:,:,Nrhs) )  ! bottom heat flux 
     242      IF( ln_trabbl  )   CALL tra_bbl       ( kstp, Nm1, Nnn, Nnn_2lev, ts(:,:,:,:,Nrhs) )  ! advective (and/or diffusive) bottom boundary layer scheme 
     243      IF( ln_tradmp  )   CALL tra_dmp       ( kstp, Nm1, Nnn_2lev, ts(:,:,:,:,Nrhs) )  ! internal damping trends 
    244244      IF( ln_bdy     )   CALL bdy_tra_dmp   ( kstp )  ! bdy damping trends 
    245245#if defined key_agrif 
     
    247247               &         CALL Agrif_Sponge_tra        ! tracers sponge 
    248248#endif 
    249                          CALL tra_adv       ( kstp, Nm1, Nnn, Np1, ts(:,:,:,:,Nrhs) )  ! horizontal & vertical advection 
     249                         CALL tra_adv       ( kstp, Nm1, Nnn, Np1, Nnn_2lev, ts(:,:,:,:,Nrhs) )  ! horizontal & vertical advection 
    250250      IF( ln_zdfosm  )   CALL tra_osm       ( kstp )  ! OSMOSIS non-local tracer fluxes 
    251251      IF( lrst_oce .AND. ln_zdfosm ) & 
    252252           &             CALL osm_rst( kstp, 'WRITE' )! write OSMOSIS outputs + wn (so must do here) to restarts 
    253                          CALL tra_ldf       ( kstp )  ! lateral mixing 
     253                         CALL tra_ldf       ( kstp, Nm1, Nnn, Nnn_2lev, ts(:,:,:,:,Nrhs) )  ! lateral mixing 
    254254 
    255255!!gm : why CALL to dia_ptr has been moved here??? (use trends info?) 
     
    343343      vb => vv(:,:,:,Nm1); vn => vv(:,:,:,Nnn); va => vv(:,:,:,Np1) 
    344344      wn => ww(:,:,:) 
     345      hdivn => hdiv(:,:,:) 
     346      rab_b  => r_ab  (:,:,:,:,Nm1_2lev); rab_n  => r_ab (:,:,:,:,Nnn_2lev) 
     347      rn2b   => r_n2 (:,:,:,Nm1_2lev)   ; rn2    => r_n2 (:,:,:,Nnn_2lev) 
    345348 
    346349      tsb => ts(:,:,:,:,Nm1); tsn => ts(:,:,:,:,Nnn); tsa => ts(:,:,:,:,Np1) 
     
    360363      gde3w_n => gde3w(:,:,:) 
    361364 
    362       ht_n => ht(:,:) 
    363       hu_b => hu(:,:,Nm1); hu_n => hu(:,:,Nnn); hu_a => hu(:,:,Np1) 
    364       hv_b => hv(:,:,Nm1); hv_n => hv(:,:,Nnn); hv_a => hv(:,:,Np1) 
    365       r1_hu_b => r1_hu(:,:,Nm1); r1_hu_n => r1_hu(:,:,Nnn); r1_hu_a => r1_hu(:,:,Np1) 
    366       r1_hv_b => r1_hv(:,:,Nm1); r1_hv_n => r1_hv(:,:,Nnn); r1_hv_a => r1_hv(:,:,Np1) 
    367  
    368365   END SUBROUTINE update_pointers 
    369366 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcadv.F90

    r10802 r10806  
    6868CONTAINS 
    6969 
    70    SUBROUTINE trc_adv( kt, ktlev1, ktlev2, ktlev3 ) 
     70   SUBROUTINE trc_adv( kt, ktlev1, ktlev2, ktlev3, kt2lev ) 
    7171      !!---------------------------------------------------------------------- 
    7272      !!                  ***  ROUTINE trc_adv  *** 
     
    7878      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    7979      INTEGER, INTENT(in) ::   ktlev1, ktlev2, ktlev3   ! time level indices for source terms 
     80      INTEGER, INTENT(in) ::   kt2lev                   ! time level index for 2-time-level source terms 
    8081      ! 
    8182      INTEGER ::   jk   ! dummy loop index 
     
    128129         CALL tra_adv_fct( kt, nittrc000, ktlev1, ktlev2, ktlev3, 'TRC', r2dttrc, zun, zvn, zwn, trb, trn, tra, jptra, nn_fct_h, nn_fct_v ) 
    129130      CASE ( np_MUS )                                 ! MUSCL 
    130          CALL tra_adv_mus( kt, nittrc000, ktlev2, 'TRC', r2dttrc, zun, zvn, zwn, trb,      tra, jptra        , ln_mus_ups )  
     131         CALL tra_adv_mus( kt, nittrc000, ktlev2, kt2lev, 'TRC', r2dttrc, zun, zvn, zwn, trb,      tra, jptra        , ln_mus_ups )  
    131132      CASE ( np_UBS )                                 ! UBS 
    132133         CALL tra_adv_ubs( kt, nittrc000, ktlev2, 'TRC', r2dttrc, zun, zvn, zwn, trb, trn, tra, jptra          , nn_ubs_v ) 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trctrp.F90

    r10802 r10806  
    6464         IF( ln_trcdmp )        CALL trc_dmp    ( kt )      ! internal damping trends 
    6565         IF( ln_bdy )           CALL trc_bdy_dmp( kt )      ! BDY damping trends 
    66                                 CALL trc_adv    ( kt, Nm1, Nnn, Np1 )      ! horizontal & vertical advection  
     66                                CALL trc_adv    ( kt, Nm1, Nnn, Np1, Nnn_2lev )      ! horizontal & vertical advection  
    6767         !                                                         ! Partial top/bottom cell: GRADh( trb )   
    6868         IF( ln_zps ) THEN 
Note: See TracChangeset for help on using the changeset viewer.