New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 4246 – NEMO

Changeset 4246


Ignore:
Timestamp:
2013-11-19T12:46:44+01:00 (10 years ago)
Author:
acc
Message:

Branch 2013/dev_r3858_NOC_ZTC, #863. Correction to divergence calculation ordering and return to single wzv routine called twice.

Location:
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r4241 r4246  
    901901      ENDIF 
    902902      ! 
     903      DO jk = 1, jpkm1 
     904         ! Correct velocities: 
     905         un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) )*umask(:,:,jk) 
     906         vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) )*vmask(:,:,jk) 
     907         ! 
     908      END DO 
     909      ! 
    903910      !                                   !* write time-spliting arrays in the restart 
    904911      IF(lrst_oce .AND.ln_bt_fw)   CALL ts_rst( kt, 'WRITE' ) 
  • branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

    r4178 r4246  
    5050 
    5151   PUBLIC   ssh_nxt    ! called by step.F90 
    52    PUBLIC   wzv_1      ! called by step.F90 
    53    PUBLIC   wzv_2      ! called by step.F90 
     52   PUBLIC   wzv        ! called by step.F90 
    5453   PUBLIC   ssh_swp    ! called by step.F90 
    5554 
     
    152151 
    153152    
    154    SUBROUTINE wzv_1( kt ) 
    155       !!---------------------------------------------------------------------- 
    156       !!                ***  ROUTINE wzv_1  *** 
     153   SUBROUTINE wzv( kt ) 
     154      !!---------------------------------------------------------------------- 
     155      !!                ***  ROUTINE wzv  *** 
    157156      !!                    
    158157      !! ** Purpose :   compute the now vertical velocity 
     
    169168      ! 
    170169      INTEGER, INTENT(in) ::   kt           ! time step 
    171       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhdiv 
     170      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d 
     171      REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv 
    172172      ! 
    173173      INTEGER             ::   ji, jj, jk   ! dummy loop indices 
     
    175175      !!---------------------------------------------------------------------- 
    176176       
    177       IF( nn_timing == 1 )  CALL timing_start('wzv_1') 
    178       ! 
    179       CALL wrk_alloc( jpi, jpj, jpk, zhdiv )  
     177      IF( nn_timing == 1 )  CALL timing_start('wzv') 
    180178      ! 
    181179      IF( kt == nit000 ) THEN 
    182180         ! 
    183181         IF(lwp) WRITE(numout,*) 
    184          IF(lwp) WRITE(numout,*) 'wzv_1 : now vertical velocity ' 
     182         IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity ' 
    185183         IF(lwp) WRITE(numout,*) '~~~~~ ' 
    186184         ! 
     
    195193      ! 
    196194      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases 
     195         CALL wrk_alloc( jpi, jpj, jpk, zhdiv )  
     196         ! 
    197197         DO jk = 1, jpkm1 
    198198            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t) 
     
    213213         END DO 
    214214         !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0 
     215         CALL wrk_dealloc( jpi, jpj, jpk, zhdiv )  
    215216      ELSE   ! z_star and linear free surface cases 
    216217         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence 
     
    224225         wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:) 
    225226#endif 
    226  
    227       ! 
    228       CALL wrk_dealloc( jpi, jpj, jpk, zhdiv )  
    229       ! 
    230       IF( nn_timing == 1 )  CALL timing_stop('wzv_1') 
    231  
    232  
    233    END SUBROUTINE wzv_1 
    234  
    235    SUBROUTINE wzv_2( kt ) 
    236       !!---------------------------------------------------------------------- 
    237       !!                ***  ROUTINE wzv_2  *** 
    238       !!                    
    239       !! ** Purpose :   compute the now vertical velocity 
    240       !! 
    241       !! ** Method  : - Using the incompressibility hypothesis, the vertical  
    242       !!      velocity is computed by integrating the horizontal divergence   
    243       !!      from the bottom to the surface minus the scale factor evolution. 
    244       !!        The boundary conditions are w=0 at the bottom (no flux) and. 
    245       !! 
    246       !! ** action  :   wn      : now vertical velocity 
    247       !! 
    248       !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
    249       !!---------------------------------------------------------------------- 
    250       ! 
    251       INTEGER, INTENT(in) ::   kt           ! time step 
    252       REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d 
    253       REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv 
    254       ! 
    255       INTEGER             ::   ji, jj, jk   ! dummy loop indices 
    256       REAL(wp)            ::   z1_2dt       ! local scalars 
    257       !!---------------------------------------------------------------------- 
    258        
    259       IF( nn_timing == 1 )  CALL timing_start('wzv_2') 
    260       ! 
    261       CALL wrk_alloc( jpi, jpj, jpk, z3d, zhdiv )  
    262       ! 
    263       IF( kt == nit000 ) THEN 
    264          ! 
    265          IF(lwp) WRITE(numout,*) 
    266          IF(lwp) WRITE(numout,*) 'wzv_2 : now vertical velocity ' 
    267          IF(lwp) WRITE(numout,*) '~~~~~ ' 
    268          ! 
    269       ENDIF 
    270       !                                           !------------------------------! 
    271       !                                           !     Now Vertical Velocity    ! 
    272       !                                           !------------------------------! 
    273       z1_2dt = 1. / ( 2. * rdt )                         ! set time step size (Euler/Leapfrog) 
    274       IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1. / rdt 
    275  
    276       IF(  lk_dynspg_ts ) THEN 
    277          ! At this stage:  
    278          ! 1) vertical scale factors have been updated.  
    279          ! 2) Time averaged barotropic velocities around step kt+1 (ua_b, va_b) as well as  
    280          !   "advective" barotropic velocities (un_adv, vn_adv) at step kt in agreement  
    281          !   with continuity equation are available. 
    282          ! 3) 3d velocities (ua, va) hold ua_b, va_b issued from time splitting as barotropic components. 
    283          ! Below => Update "now" velocities, divergence, then vertical velocity 
    284          ! NB: hdivn is NOT updated such that hdivb is not updated also. This means that horizontal 
    285          !     momentum diffusion is still performed on Instantaneous barotropic arrays. I experencied 
    286          !     some issues with UBS with the current method. See "divup" comments in development branch to  
    287          !     update the divergence fully if necessary (2013/dev_r3867_MERCATOR1_DYN). 
    288          ! 
    289          DO jk = 1, jpkm1 
    290             ! Correct velocities:                             
    291             un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) )*umask(:,:,jk) 
    292             vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) )*vmask(:,:,jk) 
    293             ! 
    294             ! Compute divergence: 
    295             DO jj = 2, jpjm1 
    296                DO ji = fs_2, fs_jpim1   ! vector opt. 
    297                   z3d(ji,jj,jk) =   & 
    298                      (  e2u(ji,jj)*fse3u(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj)*fse3u(ji-1,jj,jk) * un(ji-1,jj,jk)       & 
    299                       + e1v(ji,jj)*fse3v(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji,jj-1)*fse3v(ji,jj-1,jk) * vn(ji,jj-1,jk)  )    & 
    300                      / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    301                END DO 
    302             END DO 
    303          END DO 
    304       ! 
    305          IF( ln_rnf      )   CALL sbc_rnf_div( z3d )      ! runoffs (update divergence) 
    306       ELSE 
    307          z3d(:,:,:) = hdivn(:,:,:) 
    308       ENDIF 
    309       ! 
    310       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases 
    311          DO jk = 1, jpkm1 
    312             ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t) 
    313             ! - ML - note: computation allready done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way) 
    314             DO jj = 2, jpjm1 
    315                DO ji = fs_2, fs_jpim1   ! vector opt. 
    316                   zhdiv(ji,jj,jk) = r1_e12t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) ) 
    317                END DO 
    318             END DO 
    319          END DO 
    320          CALL lbc_lnk(zhdiv, 'T', 1.)  ! - ML - Perhaps not necessary: not used for horizontal "connexions" 
    321          !                             ! Is it problematic to have a wrong vertical velocity in boundary cells? 
    322          !                             ! Same question holds for hdivn. Perhaps just for security 
    323          DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence 
    324             ! computation of w 
    325             wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * z3d(:,:,jk) + zhdiv(:,:,jk)                    & 
    326                &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk) 
    327          END DO 
    328          !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0 
    329       ELSE   ! z_star and linear free surface cases 
    330          DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence 
    331             ! computation of w 
    332             wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * z3d(:,:,jk)                                    & 
    333                &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk) 
    334          END DO 
    335       ENDIF 
    336  
    337 #if defined key_bdy 
    338          wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:) 
    339 #endif 
    340  
     227      ! 
    341228      !                                           !------------------------------! 
    342229      !                                           !           outputs            ! 
     
    345232      IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value 
    346233         CALL wrk_alloc( jpi, jpj, z2d )  
     234         CALL wrk_alloc( jpi, jpj, jpk, z3d )  
    347235         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport. 
    348236         z2d(:,:) = rau0 * e12t(:,:) 
     
    353241         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) 
    354242         CALL wrk_dealloc( jpi, jpj, z2d  )  
    355       ENDIF 
    356       ! 
    357       CALL wrk_dealloc( jpi, jpj, jpk, z3d, zhdiv )  
    358       ! 
    359       IF( nn_timing == 1 )  CALL timing_stop('wzv_2') 
    360  
    361  
    362    END SUBROUTINE wzv_2 
     243         CALL wrk_dealloc( jpi, jpj, jpk, z3d )  
     244      ENDIF 
     245      ! 
     246      IF( nn_timing == 1 )  CALL timing_stop('wzv') 
     247 
     248 
     249   END SUBROUTINE wzv 
    363250 
    364251   SUBROUTINE ssh_swp( kt ) 
  • branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/step.F90

    r4178 r4246  
    9898      !  Ocean dynamics : hdiv, rot, ssh, e3, wn 
    9999      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    100                          CALL ssh_nxt       ( kstp )  ! after ssh 
     100                         CALL ssh_nxt       ( kstp )  ! after ssh (includes call to div_cur) 
    101101      IF( lk_dynspg_ts ) THEN 
    102                                   CALL wzv_1         ( kstp )  ! now cross-level velocity  
     102                                  CALL wzv           ( kstp )  ! now cross-level velocity  
    103103          ! In case the time splitting case, update almost all momentum trends here: 
    104104          ! Note that the computation of vertical velocity above, hence "after" sea level 
     
    124124                                  CALL dyn_hpg( kstp )         ! horizontal gradient of Hydrostatic pressure 
    125125                                  CALL dyn_spg( kstp, indic )  ! surface pressure gradient 
     126                                  CALL div_cur( kstp )         ! Horizontal divergence & Relative vorticity (2nd call in time-split case) 
    126127 
    127128                                  ua_sv(:,:,:) = ua(:,:,:)     ! save next velocities (not trends !) 
     
    129130      ENDIF 
    130131      IF( lk_vvl     )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors 
    131                          CALL wzv_2         ( kstp )  ! now cross-level velocity (original) 
     132                         CALL wzv           ( kstp )  ! now cross-level velocity (original) 
    132133 
    133134      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
  • branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

    r3865 r4246  
    1111   USE ldftra_oce       ! ocean tracer   - trends 
    1212   USE ldfdyn_oce       ! ocean dynamics - trends 
     13   USE divcur           ! hor. divergence and curl      (div & cur routines) 
    1314   USE in_out_manager   ! I/O manager 
    1415   USE iom              ! 
Note: See TracChangeset for help on using the changeset viewer.