Changeset 10919


Ignore:
Timestamp:
2019-04-30T17:29:23+02:00 (18 months ago)
Author:
davestorkey
Message:

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps :

  1. Introduce ssh, uu_b, vv_b variables with time dimension (and pointers).
  2. Convert dynspg modules (but no change to subcycling in dynspg_ts.F90).
Location:
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE
Files:
5 edited

Legend:

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

    r10068 r10919  
    5353CONTAINS 
    5454 
    55    SUBROUTINE dyn_spg( kt ) 
     55   SUBROUTINE dyn_spg( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa ) 
    5656      !!---------------------------------------------------------------------- 
    5757      !!                  ***  ROUTINE dyn_spg  *** 
     
    7171      !!             period is used to prevent the divergence of odd and even time step. 
    7272      !!---------------------------------------------------------------------- 
    73       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     73      INTEGER                             , INTENT( in )  ::  kt                  ! ocean time-step index 
     74      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs, Kaa ! ocean time level indices 
     75      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv            ! ocean velocities and RHS of momentum equation 
     76      REAL(wp), DIMENSION(jpi,jpj,jpt)    , INTENT(inout) ::  pssh, puu_b, pvv_b  ! SSH and barotropic velocities at main time levels 
    7477      ! 
    7578      INTEGER  ::   ji, jj, jk                   ! dummy loop indices 
     
    8386      IF( l_trddyn )   THEN                      ! temporary save of ta and sa trends 
    8487         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )  
    85          ztrdu(:,:,:) = ua(:,:,:) 
    86          ztrdv(:,:,:) = va(:,:,:) 
     88         ztrdu(:,:,:) = puu(:,:,:,Krhs) 
     89         ztrdv(:,:,:) = pvv(:,:,:,Krhs) 
    8790      ENDIF 
    8891      ! 
     
    126129               DO jj = 2, jpjm1                    ! add scalar approximation for load potential 
    127130                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    128                      spgu(ji,jj) = spgu(ji,jj) + zld * ( sshn(ji+1,jj) - sshn(ji,jj) ) * r1_e1u(ji,jj) 
    129                      spgv(ji,jj) = spgv(ji,jj) + zld * ( sshn(ji,jj+1) - sshn(ji,jj) ) * r1_e2v(ji,jj) 
     131                     spgu(ji,jj) = spgu(ji,jj) + zld * ( pssh(ji+1,jj,Kmm) - pssh(ji,jj,Kmm) ) * r1_e1u(ji,jj) 
     132                     spgv(ji,jj) = spgv(ji,jj) + zld * ( pssh(ji,jj+1,Kmm) - pssh(ji,jj,Kmm) ) * r1_e2v(ji,jj) 
    130133                  END DO  
    131134               END DO 
     
    150153            DO jj = 2, jpjm1 
    151154               DO ji = fs_2, fs_jpim1   ! vector opt. 
    152                   ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj) 
    153                   va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj) 
     155                  puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + spgu(ji,jj) 
     156                  pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + spgv(ji,jj) 
    154157               END DO 
    155158            END DO 
     
    161164      ! 
    162165      SELECT CASE ( nspg )                   !== surface pressure gradient computed and add to the general trend ==! 
    163       CASE ( np_EXP )   ;   CALL dyn_spg_exp( kt )              ! explicit 
    164       CASE ( np_TS  )   ;   CALL dyn_spg_ts ( kt )              ! time-splitting 
     166      CASE ( np_EXP )   ;   CALL dyn_spg_exp( kt,      Kmm,       puu, pvv, Krhs )                    ! explicit 
     167      CASE ( np_TS  )   ;   CALL dyn_spg_ts ( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa ) ! time-splitting 
    165168      END SELECT 
    166169      !                     
    167170      IF( l_trddyn )   THEN                  ! save the surface pressure gradient trends for further diagnostics 
    168          ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    169          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     171         ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:) 
     172         ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:) 
    170173         CALL trd_dyn( ztrdu, ztrdv, jpdyn_spg, kt ) 
    171174         DEALLOCATE( ztrdu , ztrdv )  
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynspg_exp.F90

    r10068 r10919  
    3838CONTAINS 
    3939 
    40    SUBROUTINE dyn_spg_exp( kt ) 
     40   SUBROUTINE dyn_spg_exp( kt, Kmm, puu, pvv, Krhs ) 
    4141      !!---------------------------------------------------------------------- 
    4242      !!                  ***  routine dyn_spg_exp  *** 
     
    4848      !! ** Method  :   Explicit free surface formulation. Add to the general 
    4949      !!              momentum trend the surface pressure gradient : 
    50       !!                      (ua,va) = (ua,va) + (spgu,spgv) 
    51       !!              where spgu = -1/rau0 d/dx(ps) = -g/e1u di( sshn ) 
    52       !!                    spgv = -1/rau0 d/dy(ps) = -g/e2v dj( sshn ) 
     50      !!                      (uu(rhs),vv(rhs)) = (uu(rhs),vv(rhs)) + (spgu,spgv) 
     51      !!              where spgu = -1/rau0 d/dx(ps) = -g/e1u di( ssh(now) ) 
     52      !!                    spgv = -1/rau0 d/dy(ps) = -g/e2v dj( ssh(now) ) 
    5353      !! 
    54       !! ** Action :   (ua,va)   trend of horizontal velocity increased by  
     54      !! ** Action :   (puu(:,:,:,Krhs),pvv(:,:,:,Krhs))   trend of horizontal velocity increased by  
    5555      !!                         the surf. pressure gradient trend 
    5656      !!--------------------------------------------------------------------- 
    57       INTEGER, INTENT(in)  ::   kt   ! ocean time-step index 
     57      INTEGER                             , INTENT( in )  ::  kt        ! ocean time-step index 
     58      INTEGER                             , INTENT( in )  ::  Kmm, Krhs ! ocean time level indices 
     59      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv  ! ocean velocities and RHS of momentum equation 
    5860      !! 
    5961      INTEGER ::   ji, jj, jk   ! dummy loop indices 
     
    7476         DO jj = 2, jpjm1                    ! now surface pressure gradient 
    7577            DO ji = fs_2, fs_jpim1   ! vector opt. 
    76                spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) * r1_e1u(ji,jj) 
    77                spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) * r1_e2v(ji,jj) 
     78               spgu(ji,jj) = - grav * ( ssh(ji+1,jj,Kmm) - ssh(ji,jj,Kmm) ) * r1_e1u(ji,jj) 
     79               spgv(ji,jj) = - grav * ( ssh(ji,jj+1,Kmm) - ssh(ji,jj,Kmm) ) * r1_e2v(ji,jj) 
    7880            END DO  
    7981         END DO 
     
    8284            DO jj = 2, jpjm1 
    8385               DO ji = fs_2, fs_jpim1   ! vector opt. 
    84                   ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj) 
    85                   va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj) 
     86                  puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + spgu(ji,jj) 
     87                  pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + spgv(ji,jj) 
    8688               END DO 
    8789            END DO 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynspg_ts.F90

    r10481 r10919  
    11MODULE dynspg_ts 
    22 
    3    !! Includes ROMS wd scheme with diagnostic outputs ; un and ua updates are commented out !  
     3   !! Includes ROMS wd scheme with diagnostic outputs ; puu(:,:,:,Kmm) and puu(:,:,:,Krhs) updates are commented out !  
    44 
    55   !!====================================================================== 
     
    119119 
    120120 
    121    SUBROUTINE dyn_spg_ts( kt ) 
     121   SUBROUTINE dyn_spg_ts( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa ) 
    122122      !!---------------------------------------------------------------------- 
    123123      !! 
     
    134134      !! 
    135135      !! ** Action : 
    136       !!      -Update the filtered free surface at step "n+1"      : ssha 
    137       !!      -Update filtered barotropic velocities at step "n+1" : ua_b, va_b 
     136      !!      -Update the filtered free surface at step "n+1"      : pssh(:,:,Kaa) 
     137      !!      -Update filtered barotropic velocities at step "n+1" : puu_b(:,:,:,Kaa), vv_b(:,:,:,Kaa) 
    138138      !!      -Compute barotropic advective fluxes at step "n"     : un_adv, vn_adv 
    139139      !!      These are used to advect tracers and are compliant with discrete 
    140140      !!      continuity equation taken at the baroclinic time steps. This  
    141141      !!      ensures tracers conservation. 
    142       !!      - (ua, va) momentum trend updated with barotropic component. 
     142      !!      - (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) momentum trend updated with barotropic component. 
    143143      !! 
    144144      !! References : Shchepetkin and McWilliams, Ocean Modelling, 2005.  
    145145      !!--------------------------------------------------------------------- 
    146       INTEGER, INTENT(in)  ::   kt   ! ocean time-step index 
     146      INTEGER                             , INTENT( in )  ::  kt                  ! ocean time-step index 
     147      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs, Kaa ! ocean time level indices 
     148      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv            ! ocean velocities and RHS of momentum equation 
     149      REAL(wp), DIMENSION(jpi,jpj,jpt)    , INTENT(inout) ::  pssh, puu_b, pvv_b  ! SSH and barotropic velocities at main time levels 
    147150      ! 
    148151      INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices 
     
    328331            DO jk = 1, jpkm1 
    329332               DO jj = 1, jpjm1 
    330                   zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 
     333                  zhf(:,jj) = zhf(:,jj) + e3f(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 
    331334               END DO 
    332335            END DO 
     
    360363      ! 
    361364      DO jk = 1, jpkm1 
    362          zu_frc(:,:) = zu_frc(:,:) + e3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 
    363          zv_frc(:,:) = zv_frc(:,:) + e3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)          
     365         zu_frc(:,:) = zu_frc(:,:) + e3u(:,:,jk,Kmm) * puu(:,:,jk,Krhs) * umask(:,:,jk) 
     366         zv_frc(:,:) = zv_frc(:,:) + e3v(:,:,jk,Kmm) * pvv(:,:,jk,Krhs) * vmask(:,:,jk)          
    364367      END DO 
    365368      ! 
     
    372375         DO jj = 2, jpjm1 
    373376            DO ji = fs_2, fs_jpim1   ! vector opt. 
    374                ua(ji,jj,jk) = ua(ji,jj,jk) - zu_frc(ji,jj) * umask(ji,jj,jk) 
    375                va(ji,jj,jk) = va(ji,jj,jk) - zv_frc(ji,jj) * vmask(ji,jj,jk) 
     377               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - zu_frc(ji,jj) * umask(ji,jj,jk) 
     378               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - zv_frc(ji,jj) * vmask(ji,jj,jk) 
    376379            END DO 
    377380         END DO 
     
    385388      !                                   ! -------------------------------------------------------- 
    386389      ! 
    387       zwx(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:)        ! now fluxes  
    388       zwy(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:) 
     390      zwx(:,:) = puu_b(:,:,Kmm) * hu_n(:,:) * e2u(:,:)        ! now fluxes  
     391      zwy(:,:) = pvv_b(:,:,Kmm) * hv_n(:,:) * e1v(:,:) 
    389392      ! 
    390393      SELECT CASE( nvor_scheme ) 
     
    393396            DO ji = 2, jpim1   ! vector opt. 
    394397               zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * r1_hu_n(ji,jj)                    & 
    395                   &               * (  e1e2t(ji+1,jj)*ht_n(ji+1,jj)*ff_t(ji+1,jj) * ( vn_b(ji+1,jj) + vn_b(ji+1,jj-1) )   & 
    396                   &                  + e1e2t(ji  ,jj)*ht_n(ji  ,jj)*ff_t(ji  ,jj) * ( vn_b(ji  ,jj) + vn_b(ji  ,jj-1) )   ) 
     398                  &               * (  e1e2t(ji+1,jj)*ht_n(ji+1,jj)*ff_t(ji+1,jj) * ( pvv_b(ji+1,jj,Kmm) + pvv_b(ji+1,jj-1,Kmm) )   & 
     399                  &                  + e1e2t(ji  ,jj)*ht_n(ji  ,jj)*ff_t(ji  ,jj) * ( pvv_b(ji  ,jj,Kmm) + pvv_b(ji  ,jj-1,Kmm) )   ) 
    397400                  ! 
    398401               zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * r1_hv_n(ji,jj)                    & 
    399                   &               * (  e1e2t(ji,jj+1)*ht_n(ji,jj+1)*ff_t(ji,jj+1) * ( un_b(ji,jj+1) + un_b(ji-1,jj+1) )   &  
    400                   &                  + e1e2t(ji,jj  )*ht_n(ji,jj  )*ff_t(ji,jj  ) * ( un_b(ji,jj  ) + un_b(ji-1,jj  ) )   )  
     402                  &               * (  e1e2t(ji,jj+1)*ht_n(ji,jj+1)*ff_t(ji,jj+1) * ( puu_b(ji,jj+1,Kmm) + puu_b(ji-1,jj+1,Kmm) )   &  
     403                  &                  + e1e2t(ji,jj  )*ht_n(ji,jj  )*ff_t(ji,jj  ) * ( puu_b(ji,jj  ,Kmm) + puu_b(ji-1,jj  ,Kmm) )   )  
    401404            END DO   
    402405         END DO   
     
    449452            DO jj = 2, jpjm1 
    450453               DO ji = 2, jpim1  
    451                   ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji+1,jj) ) >                & 
     454                  ll_tmp1 = MIN(  pssh(ji,jj,Kmm)               ,  pssh(ji+1,jj,Kmm) ) >                & 
    452455                     &      MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            & 
    453                      &      MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji+1,jj) + ht_0(ji+1,jj) )  & 
     456                     &      MAX(  pssh(ji,jj,Kmm) + ht_0(ji,jj) ,  pssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) )  & 
    454457                     &                                                         > rn_wdmin1 + rn_wdmin2 
    455                   ll_tmp2 = ( ABS( sshn(ji+1,jj)            -  sshn(ji  ,jj))  > 1.E-12 ).AND.( & 
    456                      &      MAX(   sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                & 
     458                  ll_tmp2 = ( ABS( pssh(ji+1,jj,Kmm)            -  pssh(ji  ,jj,Kmm))  > 1.E-12 ).AND.( & 
     459                     &      MAX(   pssh(ji,jj,Kmm)              ,  pssh(ji+1,jj,Kmm) ) >                & 
    457460                     &      MAX(  -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
    458461                  IF(ll_tmp1) THEN 
    459462                     zcpx(ji,jj) = 1.0_wp 
    460463                  ELSEIF(ll_tmp2) THEN 
    461                      ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
    462                      zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 
    463                                  &    / (sshn(ji+1,jj) - sshn(ji  ,jj)) ) 
     464                     ! no worries about  pssh(ji+1,jj,Kmm) -  pssh(ji  ,jj,Kmm) = 0, it won't happen ! here 
     465                     zcpx(ji,jj) = ABS( (pssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - pssh(ji,jj,Kmm) - ht_0(ji,jj)) & 
     466                                 &    / (pssh(ji+1,jj,Kmm) - pssh(ji  ,jj,Kmm)) ) 
    464467                     zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 
    465468                  ELSE 
     
    467470                  ENDIF 
    468471                  ! 
    469                   ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji,jj+1) ) >                & 
     472                  ll_tmp1 = MIN(  pssh(ji,jj,Kmm)               ,  pssh(ji,jj+1,Kmm) ) >                & 
    470473                     &      MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) .AND.            & 
    471                      &      MAX(  sshn(ji,jj) + ht_0(ji,jj) ,  sshn(ji,jj+1) + ht_0(ji,jj+1) )  & 
     474                     &      MAX(  pssh(ji,jj,Kmm) + ht_0(ji,jj) ,  pssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) )  & 
    472475                     &                                                       > rn_wdmin1 + rn_wdmin2 
    473                   ll_tmp2 = ( ABS( sshn(ji,jj)              -  sshn(ji,jj+1))  > 1.E-12 ).AND.( & 
    474                      &      MAX(   sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                & 
     476                  ll_tmp2 = ( ABS( pssh(ji,jj,Kmm)              -  pssh(ji,jj+1,Kmm))  > 1.E-12 ).AND.( & 
     477                     &      MAX(   pssh(ji,jj,Kmm)              ,  pssh(ji,jj+1,Kmm) ) >                & 
    475478                     &      MAX(  -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
    476479   
     
    478481                     zcpy(ji,jj) = 1.0_wp 
    479482                  ELSE IF(ll_tmp2) THEN 
    480                      ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here 
    481                      zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 
    482                         &             / (sshn(ji,jj+1) - sshn(ji,jj  )) ) 
     483                     ! no worries about  pssh(ji,jj+1,Kmm) -  pssh(ji,jj  ,Kmm) = 0, it won't happen ! here 
     484                     zcpy(ji,jj) = ABS( (pssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - pssh(ji,jj,Kmm) - ht_0(ji,jj)) & 
     485                        &             / (pssh(ji,jj+1,Kmm) - pssh(ji,jj  ,Kmm)) ) 
    483486                     zcpy(ji,jj) = MAX(  0._wp , MIN( zcpy(ji,jj) , 1.0_wp )  ) 
    484487                  ELSE 
     
    490493            DO jj = 2, jpjm1 
    491494               DO ji = 2, jpim1 
    492                   zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   & 
     495                  zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( pssh(ji+1,jj  ,Kmm) - pssh(ji  ,jj ,Kmm) )   & 
    493496                     &                          * r1_e1u(ji,jj) * zcpx(ji,jj)  * wdrampu(ji,jj)  !jth 
    494                   zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) )   & 
     497                  zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( pssh(ji  ,jj+1,Kmm) - pssh(ji  ,jj ,Kmm) )   & 
    495498                     &                          * r1_e2v(ji,jj) * zcpy(ji,jj)  * wdrampv(ji,jj)  !jth 
    496499               END DO 
     
    501504            DO jj = 2, jpjm1 
    502505               DO ji = fs_2, fs_jpim1   ! vector opt. 
    503                   zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) * r1_e1u(ji,jj) 
    504                   zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) * r1_e2v(ji,jj)  
     506                  zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  pssh(ji+1,jj  ,Kmm) - pssh(ji  ,jj  ,Kmm)  ) * r1_e1u(ji,jj) 
     507                  zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  pssh(ji  ,jj+1,Kmm) - pssh(ji  ,jj  ,Kmm)  ) * r1_e2v(ji,jj)  
    505508               END DO 
    506509            END DO 
     
    522525               ikbu = mbku(ji,jj)        
    523526               ikbv = mbkv(ji,jj)     
    524                zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities 
    525                zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj) 
     527               zwx(ji,jj) = puu(ji,jj,ikbu,Kmm) - puu_b(ji,jj,Kmm) ! NOW bottom baroclinic velocities 
     528               zwy(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm) 
    526529            END DO 
    527530         END DO 
     
    531534               ikbu = mbku(ji,jj)        
    532535               ikbv = mbkv(ji,jj)     
    533                zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities 
    534                zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj) 
     536               zwx(ji,jj) = puu(ji,jj,ikbu,Kbb) - puu_b(ji,jj,Kbb) ! BEFORE bottom baroclinic velocities 
     537               zwy(ji,jj) = pvv(ji,jj,ikbv,Kbb) - pvv_b(ji,jj,Kbb) 
    535538            END DO 
    536539         END DO 
     
    563566                  iktu = miku(ji,jj) 
    564567                  iktv = mikv(ji,jj) 
    565                   zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities 
    566                   zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj) 
     568                  zwx(ji,jj) = puu(ji,jj,iktu,Kmm) - puu_b(ji,jj,Kmm) ! NOW top baroclinic velocities 
     569                  zwy(ji,jj) = pvv(ji,jj,iktv,Kmm) - pvv_b(ji,jj,Kmm) 
    567570               END DO 
    568571            END DO 
     
    572575                  iktu = miku(ji,jj) 
    573576                  iktv = mikv(ji,jj) 
    574                   zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities 
    575                   zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj) 
     577                  zwx(ji,jj) = puu(ji,jj,iktu,Kbb) - puu_b(ji,jj,Kbb) ! BEFORE top baroclinic velocities 
     578                  zwy(ji,jj) = pvv(ji,jj,iktv,Kbb) - pvv_b(ji,jj,Kbb) 
    576579               END DO 
    577580            END DO 
     
    674677      ! 
    675678      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                     
    676          sshn_e(:,:) =    sshn(:,:)             
    677          un_e  (:,:) =    un_b(:,:)             
    678          vn_e  (:,:) =    vn_b(:,:) 
     679         sshn_e(:,:) =    pssh(:,:,Kmm)             
     680         un_e  (:,:) =    puu_b(:,:,Kmm)             
     681         vn_e  (:,:) =    pvv_b(:,:,Kmm) 
    679682         ! 
    680683         hu_e  (:,:) =    hu_n(:,:)        
     
    683686         hvr_e (:,:) = r1_hv_n(:,:) 
    684687      ELSE                                ! CENTRED integration: start from BEFORE fields 
    685          sshn_e(:,:) =    sshb(:,:) 
    686          un_e  (:,:) =    ub_b(:,:)          
    687          vn_e  (:,:) =    vb_b(:,:) 
     688         sshn_e(:,:) =    pssh(:,:,Kbb) 
     689         un_e  (:,:) =    puu_b(:,:,Kbb)          
     690         vn_e  (:,:) =    pvv_b(:,:,Kbb) 
    688691         ! 
    689692         hu_e  (:,:) =    hu_b(:,:)        
     
    696699      ! 
    697700      ! Initialize sums: 
    698       ua_b  (:,:) = 0._wp       ! After barotropic velocities (or transport if flux form)           
    699       va_b  (:,:) = 0._wp 
    700       ssha  (:,:) = 0._wp       ! Sum for after averaged sea level 
     701      puu_b  (:,:,Kaa) = 0._wp       ! After barotropic velocities (or transport if flux form)           
     702      pvv_b  (:,:,Kaa) = 0._wp 
     703      pssh  (:,:,Kaa) = 0._wp       ! Sum for after averaged sea level 
    701704      un_adv(:,:) = 0._wp       ! Sum for now transport issued from ts loop 
    702705      vn_adv(:,:) = 0._wp 
     
    11901193         za1 = wgtbtp1(jn)                                     
    11911194         IF( ln_dynadv_vec .OR. ln_linssh ) THEN    ! Sum velocities 
    1192             ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:)  
    1193             va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:)  
     1195            puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:)  
     1196            pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:)  
    11941197         ELSE                                       ! Sum transports 
    11951198            IF ( .NOT.ln_wd_dl ) THEN   
    1196                ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:) 
    1197                va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:) 
     1199               puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:) * hu_e (:,:) 
     1200               pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:) * hv_e (:,:) 
    11981201            ELSE  
    1199                ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:) * zuwdmask(:,:) 
    1200                va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:) * zvwdmask(:,:) 
     1202               puu_b  (:,:,Kaa) = puu_b  (:,:,Kaa) + za1 * ua_e  (:,:) * hu_e (:,:) * zuwdmask(:,:) 
     1203               pvv_b  (:,:,Kaa) = pvv_b  (:,:,Kaa) + za1 * va_e  (:,:) * hv_e (:,:) * zvwdmask(:,:) 
    12011204            END IF  
    12021205         ENDIF 
    12031206         !                                          ! Sum sea level 
    1204          ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:) 
     1207         pssh(:,:,Kaa) = pssh(:,:,Kaa) + za1 * ssha_e(:,:) 
    12051208 
    12061209         !                                                 ! ==================== ! 
     
    12361239      IF( ln_dynadv_vec .OR. ln_linssh ) THEN 
    12371240         DO jk=1,jpkm1 
    1238             ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * r1_2dt_b 
    1239             va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * r1_2dt_b 
     1241            puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) ) * r1_2dt_b 
     1242            pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) ) * r1_2dt_b 
    12401243         END DO 
    12411244      ELSE 
    1242          ! At this stage, ssha has been corrected: compute new depths at velocity points 
     1245         ! At this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points 
    12431246         DO jj = 1, jpjm1 
    12441247            DO ji = 1, jpim1      ! NO Vector Opt. 
    12451248               zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj)  * r1_e1e2u(ji,jj) & 
    1246                   &              * ( e1e2t(ji  ,jj) * ssha(ji  ,jj)      & 
    1247                   &              +   e1e2t(ji+1,jj) * ssha(ji+1,jj) ) 
     1249                  &              * ( e1e2t(ji  ,jj) * pssh(ji  ,jj,Kaa)      & 
     1250                  &              +   e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) ) 
    12481251               zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj)  * r1_e1e2v(ji,jj) & 
    1249                   &              * ( e1e2t(ji,jj  ) * ssha(ji,jj  )      & 
    1250                   &              +   e1e2t(ji,jj+1) * ssha(ji,jj+1) ) 
     1252                  &              * ( e1e2t(ji,jj  ) * pssh(ji,jj  ,Kaa)      & 
     1253                  &              +   e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) ) 
    12511254            END DO 
    12521255         END DO 
     
    12541257         ! 
    12551258         DO jk=1,jpkm1 
    1256             ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * r1_2dt_b 
    1257             va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * r1_2dt_b 
     1259            puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu_n(:,:) * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu_b(:,:) ) * r1_2dt_b 
     1260            pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv_n(:,:) * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv_b(:,:) ) * r1_2dt_b 
    12581261         END DO 
    12591262         ! Save barotropic velocities not transport: 
    1260          ua_b(:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 
    1261          va_b(:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 
     1263         puu_b(:,:,Kaa) =  puu_b(:,:,Kaa) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 
     1264         pvv_b(:,:,Kaa) =  pvv_b(:,:,Kaa) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 
    12621265      ENDIF 
    12631266 
     
    12651268      ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases)   
    12661269      DO jk = 1, jpkm1 
    1267          un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:)*r1_hu_n(:,:) - un_b(:,:) ) * umask(:,:,jk) 
    1268          vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:)*r1_hv_n(:,:) - vn_b(:,:) ) * vmask(:,:,jk) 
     1270         puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) + un_adv(:,:)*r1_hu_n(:,:) - puu_b(:,:,Kmm) ) * umask(:,:,jk) 
     1271         pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) + vn_adv(:,:)*r1_hv_n(:,:) - pvv_b(:,:,Kmm) ) * vmask(:,:,jk) 
    12691272      END DO 
    12701273 
    12711274      IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN  
    12721275         DO jk = 1, jpkm1 
    1273             un(:,:,jk) = ( un_adv(:,:)*r1_hu_n(:,:) & 
    1274                        & + zuwdav2(:,:)*(un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:)) ) * umask(:,:,jk)  
    1275             vn(:,:,jk) = ( vn_adv(:,:)*r1_hv_n(:,:) &  
    1276                        & + zvwdav2(:,:)*(vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:)) ) * vmask(:,:,jk)   
     1276            puu(:,:,jk,Kmm) = ( un_adv(:,:)*r1_hu_n(:,:) & 
     1277                       & + zuwdav2(:,:)*(puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu_n(:,:)) ) * umask(:,:,jk)  
     1278            pvv(:,:,jk,Kmm) = ( vn_adv(:,:)*r1_hv_n(:,:) &  
     1279                       & + zvwdav2(:,:)*(pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv_n(:,:)) ) * vmask(:,:,jk)   
    12771280         END DO 
    12781281      END IF  
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/oce.F90

    r10876 r10919  
    3333   !! free surface                                      !  before  ! now    ! after  ! 
    3434   !! ------------                                      !  fields  ! fields ! fields ! 
    35    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ub_b   ,  un_b  ,  ua_b  !: Barotropic velocities at u-point [m/s] 
    36    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vb_b   ,  vn_b  ,  va_b  !: Barotropic velocities at v-point [m/s] 
    37    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshb   ,  sshn  ,  ssha  !: sea surface height at t-point [m] 
     35   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET ::   ssh, uu_b,  vv_b   !: SSH [m] and barotropic velocities [m/s] 
    3836 
    3937   !! Arrays at barotropic time step:                   ! befbefore! before !  now   ! after  ! 
     
    6967 
    7068   !! TEMPORARY POINTERS FOR DEVELOPMENT ONLY 
    71    REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:)   ::   ub   ,  un    , ua     !: i-horizontal velocity        [m/s] 
    72    REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:)   ::   vb   ,  vn    , va     !: j-horizontal velocity        [m/s] 
    73    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(:,:,:,:) ::   tsb  ,  tsn   , tsa    !: 4D T-S fields                [Celsius,psu]          
     69   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:)   ::   ub   ,  un    , ua       !: i-horizontal velocity        [m/s] 
     70   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:)   ::   vb   ,  vn    , va       !: j-horizontal velocity        [m/s] 
     71   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:)   ::           wn               !: k-vertical   velocity        [m/s] 
     72   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:)   ::           hdivn            !: horizontal divergence        [s-1] 
     73   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:,:) ::   tsb  ,  tsn   , tsa      !: 4D T-S fields                [Celsius,psu]          
     74   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:)     ::   ub_b   ,  un_b  ,  ua_b  !: Barotropic velocities at u-point [m/s] 
     75   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:)     ::   vb_b   ,  vn_b  ,  va_b  !: Barotropic velocities at v-point [m/s] 
     76   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:)     ::   sshb   ,  sshn  ,  ssha  !: sea surface height at t-point [m] 
    7677   !! TEMPORARY POINTERS FOR DEVELOPMENT ONLY 
    7778 
     
    9899         &      rhd  (jpi,jpj,jpk)      , rhop (jpi,jpj,jpk)                              , STAT=ierr(1) ) 
    99100         ! 
    100       ALLOCATE( sshb(jpi,jpj)     , sshn(jpi,jpj)   , ssha(jpi,jpj)   ,     & 
    101          &      ub_b(jpi,jpj)     , un_b(jpi,jpj)   , ua_b(jpi,jpj)   ,     & 
    102          &      vb_b(jpi,jpj)     , vn_b(jpi,jpj)   , va_b(jpi,jpj)   ,     & 
     101      ALLOCATE( ssh(jpi,jpj,jpt)  , uu_b(jpi,jpj,jpt) , vv_b(jpi,jpj,jpt) , & 
    103102         &      spgu  (jpi,jpj)   , spgv(jpi,jpj)                     ,     & 
    104103         &      gtsu(jpi,jpj,jpts), gtsv(jpi,jpj,jpts)                ,     & 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90

    r10905 r10919  
    198198      IF( ln_zdfosm  )   CALL dyn_osm( kstp,      Nnn      , uu, vv, Nrhs )  ! OSMOSIS non-local velocity fluxes ==> RHS 
    199199                         CALL dyn_hpg       ( kstp )  ! horizontal gradient of Hydrostatic pressure 
    200                          CALL dyn_spg       ( kstp )  ! surface pressure gradient 
     200                         CALL dyn_spg       ( kstp, Nbb, Nnn, Nrhs, uu, vv, ssh, uu_b, vv_b, Naa )  ! surface pressure gradient 
    201201 
    202202                                                      ! With split-explicit free surface, since now transports have been updated and ssha as well 
     
    351351      hdivn => hdiv(:,:,:) 
    352352 
     353      sshb =>  ssh(:,:,Kbb); sshn =>  ssh(:,:,Kmm); ssha =>  ssh(:,:,Kaa) 
     354      ub_b => uu_b(:,:,Kbb); un_b => uu_b(:,:,Kmm); ua_b => uu_b(:,:,Kaa) 
     355      vb_b => vv_b(:,:,Kbb); vn_b => vv_b(:,:,Kmm); va_b => vv_b(:,:,Kaa) 
     356 
    353357      tsb => ts(:,:,:,:,Kbb); tsn => ts(:,:,:,:,Kmm); tsa => ts(:,:,:,:,Kaa) 
    354358 
Note: See TracChangeset for help on using the changeset viewer.