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

Changeset 1368


Ignore:
Timestamp:
2009-04-01T14:02:17+02:00 (15 years ago)
Author:
rblod
Message:

dev_004_VVL: change time stepping for ssh, see ticket #394

Location:
branches/dev_004_VVL/NEMO/OPA_SRC
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynspg_flt.F90

    r1200 r1368  
    427427      END DO 
    428428 
    429       ! Sea surface elevation time stepping 
    430       ! ----------------------------------- 
    431       IF( lk_vvl ) THEN   ! after free surface elevation 
    432          zssha(:,:) = ssha(:,:) 
    433       ELSE 
    434          zssha(:,:) = sshb(:,:) + z2dt * ( wn(:,:,1) - emp(:,:) * zraur ) * tmask(:,:,1) 
    435       ENDIF 
    436 #if defined key_obc 
    437 # if defined key_agrif 
    438       IF ( Agrif_Root() ) THEN 
    439 # endif 
    440          zssha(:,:)=zssha(:,:)*obctmsk(:,:) 
    441          CALL lbc_lnk(zssha,'T',1.)  ! absolutly compulsory !! (jmm) 
    442 # if defined key_agrif 
    443       ENDIF 
    444 # endif 
    445 #endif 
    446  
    447       IF( neuler == 0 .AND. kt == nit000 ) THEN      ! Euler (forward) time stepping, no time filter 
    448          ! swap of arrays 
    449          sshb(:,:) = sshn (:,:) 
    450          sshn(:,:) = zssha(:,:) 
    451       ELSE                                           ! Leap-frog time stepping and time filter 
    452          ! time filter and array swap 
    453          sshb(:,:) = atfp * ( sshb(:,:) + zssha(:,:) ) + atfp1 * sshn(:,:) 
    454          sshn(:,:) = zssha(:,:) 
    455       ENDIF 
    456  
    457429      ! write filtered free surface arrays in restart file 
    458430      ! -------------------------------------------------- 

  • branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r1241 r1368  
    566566      ! Phase 3 : Update sea surface height from time averaged barotropic variables 
    567567      ! --------------------------------------------------------------------------- 
    568  
    569   
    570       ! Horizontal divergence of time averaged barotropic transports 
    571       !------------------------------------------------------------- 
    572       IF( .NOT. lk_vvl ) THEN 
    573          DO jj = 2, jpjm1 
    574             DO ji = fs_2, fs_jpim1   ! vector opt. 
    575                zhdiv(ji,jj) = ( e2u(ji,jj) * un_b(ji,jj) - e2u(ji-1,jj  ) * un_b(ji-1,jj  )     & 
    576                   &            +e1v(ji,jj) * vn_b(ji,jj) - e1v(ji  ,jj-1) * vn_b(ji  ,jj-1) )   & 
    577                   &           / ( e1t(ji,jj) * e2t(ji,jj) ) 
    578             END DO 
    579          END DO 
    580       ENDIF 
    581  
    582 #if defined key_obc && ! defined key_vvl 
    583       ! open boundaries (div must be zero behind the open boundary) 
    584       !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column 
    585       IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1)   = 0.e0    ! east 
    586       IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1)   = 0.e0    ! west 
    587       IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0    ! north 
    588       IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0    ! south 
    589 #endif 
    590  
    591 #if defined key_bdy 
    592          DO jj = 1, jpj 
    593            DO ji = 1, jpi 
    594              zhdiv(ji,jj) = zhdiv(ji,jj)*bdytmask(ji,jj) 
    595            END DO 
    596          END DO 
    597 #endif 
    598  
    599       ! sea surface height 
    600       !------------------- 
    601       IF( lk_vvl ) THEN 
    602          sshbb(:,:) = sshb(:,:) 
    603          sshb (:,:) = sshn(:,:) 
    604          sshn (:,:) = ssha(:,:) 
    605       ELSE 
    606          sshb (:,:) = sshn(:,:) 
    607          sshn(:,:) = (  sshb_b(:,:) - z2dt_b * ( zraur * emp(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1) 
    608       ENDIF 
    609  
    610       ! ... Boundary conditions on sshn 
    611       IF( .NOT. lk_obc ) CALL lbc_lnk( sshn, 'T', 1. ) 
    612  
     568!RB_vvl now done in ssh_wzv and ssh_nxt 
    613569 
    614570      ! ----------------------------------------------------------------------------- 

  • branches/dev_004_VVL/NEMO/OPA_SRC/DYN/wzvmod.F90

    r1362 r1368  
    2929   !! * Routine accessibility 
    3030   PUBLIC wzv          ! routine called by step.F90 and inidtr.F90 
     31   PUBLIC ssh_wzv 
     32   PUBLIC ssh_nxt 
    3133 
    3234   !! * Substitutions 
    … …  
    5456      !! ** action  :   wn array : the now vertical velocity 
    5557      !!---------------------------------------------------------------------- 
    56       !! * Arguments 
    57       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
    58  
    59       !! * Local declarations 
    60       INTEGER  ::           jk           ! dummy loop indices 
    61       !! Variable volume 
    62       INTEGER  ::   ji, jj               ! dummy loop indices 
    63       REAL(wp) ::   z2dt, zraur          ! temporary scalar 
    64       REAL(wp), DIMENSION (jpi,jpj) ::   zssha, zun, zvn, zhdiv 
    65 #if defined key_bdy 
    66       INTEGER  ::     jgrd, jb           ! temporary scalars 
    67 #endif 
     58      INTEGER, INTENT(in) :: kt 
     59       
     60      ! Empty routine 
     61 
     62      WRITE(*,*) 'wzv : you should not be here : error ?'   
     63 
     64   END SUBROUTINE wzv 
     65 
     66   SUBROUTINE ssh_wzv( kt )  
     67      !!---------------------------------------------------------------------- 
     68      !!                ***  ROUTINE dom_vvl_ssh  *** 
     69      !!                    
     70      !! ** Purpose :   compute the after ssh (ssha), the now vertical velocity 
     71      !!              Cand update the now vertical coordinate (lk_vvl=T). 
     72      !! 
     73      !! ** Method  : -   
     74 
     75      !!              - Using the incompressibility hypothesis, the vertical 
     76      !!      velocity is computed by integrating the horizontal divergence  
     77      !!      from the bottom to the surface minus the scale factor evolution. 
     78      !!        The boundary conditions are w=0 at the bottom (no flux) and. 
     79      !! 
     80      !! ** action  :   wn array : the now vertical velocity 
     81      !!---------------------------------------------------------------------- 
     82      INTEGER, INTENT(in) ::   kt   ! time step 
     83      !! 
     84      INTEGER  ::   ji, jj, jk      ! dummy loop indices 
     85      REAL(wp) ::   z2dt, zraur     ! temporary scalars 
     86      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv       ! 2D workspace 
    6887      !!---------------------------------------------------------------------- 
    6988 
    7089      IF( kt == nit000 ) THEN 
    7190         IF(lwp) WRITE(numout,*) 
    72          IF(lwp) WRITE(numout,*) 'wzv     : vertical velocity from continuity eq.' 
    73          IF(lwp) WRITE(numout,*) '~~~~~~~ '  
    74  
    75          ! bottom boundary condition: w=0 (set once for all) 
    76          wn(:,:,jpk) = 0.e0 
    77       ENDIF 
    78  
    79       IF( lk_vvl ) THEN                ! Variable volume 
     91         IF(lwp) WRITE(numout,*) 'ssh_wzv : vertical velocity from continuity eq. (vvl option)' 
     92         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
    8093         ! 
    81          z2dt = 2. * rdt                                       ! time step: leap-frog 
    82          IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt       ! time step: Euler if restart from rest 
    83          zraur  = 1. / rauw 
    84  
    85          ! Vertically integrated quantities 
    86          ! -------------------------------- 
    87          zun(:,:) = 0.e0 
    88          zvn(:,:) = 0.e0 
    89          ! 
    90          DO jk = 1, jpkm1             ! Vertically integrated transports (now) 
    91             zun(:,:) = zun(:,:) + fse3u(:,:,jk) * un(:,:,jk) 
    92             zvn(:,:) = zvn(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 
     94         wn(:,:,jpk) = 0.e0      ! bottom boundary condition: w=0 (set once for all) 
     95      ENDIF 
     96 
     97      ! set time step size (Euler/Leapfrog) 
     98      z2dt = 2. * rdt  
     99      IF( neuler == 0 .AND. kt == nit000 ) z2dt =rdt 
     100 
     101      zraur = 1. / rauw 
     102 
     103      !                                           !------------------------------! 
     104      !                                           !   After Sea Surface Height   ! 
     105      !                                           !------------------------------! 
     106      zhdiv(:,:) = 0.e0 
     107      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports 
     108        zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk) 
     109      END DO 
     110 
     111      !                                                ! Sea surface elevation time stepping 
     112      ssha(:,:) = (  sshb(:,:) - z2dt * ( zraur * emp(:,:)                       + zhdiv(:,:) )  ) * tmask(:,:,1) 
     113 
     114      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only) 
     115      IF( lk_vvl ) THEN                                ! (required only in key_vvl case) 
     116!RB_vvl 
     117! Compute ssh at u,v, f points and update vertical coordinate 
     118! Currently done in dom_vvl 
     119      ENDIF 
     120 
     121      !                                           !------------------------------! 
     122      !                                           !     Now Vertical Velocity    ! 
     123      !                                           !------------------------------! 
     124      !                                                ! integrate from the bottom the hor. divergence 
     125         DO jk = jpkm1, 1, -1 
     126            wn(:,:,jk) = wn(:,:,jk+1) -    fse3t_n(:,:,jk) * hdivn(:,:,jk)   & 
     127                 &                    - (  fse3t_a(:,:,jk)                   & 
     128                 &                       - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt 
    93129         END DO 
    94130 
    95          ! Horizontal divergence of barotropic transports 
    96          !-------------------------------------------------- 
    97          zhdiv(:,:) = 0.e0 
    98          DO jj = 2, jpjm1 
    99             DO ji = 2, jpim1   ! vector opt. 
    100                zhdiv(ji,jj) = (  e2u(ji  ,jj  ) * zun(ji  ,jj  )     & 
    101                   &            - e2u(ji-1,jj  ) * zun(ji-1,jj  )     & 
    102                   &            + e1v(ji  ,jj  ) * zvn(ji  ,jj  )     & 
    103                   &            - e1v(ji  ,jj-1) * zvn(ji  ,jj-1) )   & 
    104                   &           / ( e1t(ji,jj) * e2t(ji,jj) ) 
     131      !                                           !------------------------------! 
     132      !                                           !  Update Now Vertical coord.  ! 
     133      !                                           !------------------------------! 
     134      IF( lk_vvl ) THEN                           ! only in vvl case) 
     135      !                                                ! now local depth and scale factors (stored in fse3. arrays) 
     136!RB_vvl to be done 
     137 
     138 
     139      ENDIF 
     140      ! 
     141   END SUBROUTINE ssh_wzv 
     142 
     143 
     144   SUBROUTINE ssh_nxt( kt ) 
     145      !!---------------------------------------------------------------------- 
     146      !!                    ***  ROUTINE ssh_nxt  *** 
     147      !! 
     148      !! ** Purpose :   achieve the sea surface  height time stepping by  
     149      !!              applying Asselin time filter and swapping the arrays 
     150      !!              ssha  already computed in ssh_wzv   
     151      !! 
     152      !! ** Method  : - apply Asselin time fiter to now ssh and swap : 
     153      !!             sshn = ssha + atfp * ( sshb -2 sshn + ssha ) 
     154      !!             sshn = ssha 
     155      !! 
     156      !! ** action  : - sshb, sshn   : before & now sea surface height 
     157      !!                               ready for the next time step 
     158      !!---------------------------------------------------------------------- 
     159      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
     160      !! 
     161      INTEGER  ::   ji, jj, jk           ! dummy loop indices 
     162      !!---------------------------------------------------------------------- 
     163 
     164      IF( kt == nit000 ) THEN 
     165         IF(lwp) WRITE(numout,*) 
     166         IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)' 
     167         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
     168      ENDIF 
     169 
     170      ! Time filter and swap of the ssh 
     171      ! ------------------------------- 
     172 
     173!RB_vvl ssh at u, f,v point to be added 
     174 
     175     IF( neuler == 0 .AND. kt == nit000 ) THEN   ! Euler time stepping 
     176         DO jj = 1, jpj 
     177            DO ji = 1, jpi 
     178               ! before <-- now 
     179               sshb  (ji,jj) = sshn(ji,jj)  
     180               ! now <-- after 
     181               sshn  (ji,jj) = ssha  (ji,jj) 
    105182            END DO 
    106183         END DO 
    107  
    108 #if defined key_obc && ( defined key_dynspg_exp || defined key_dynspg_ts ) 
    109          ! open boundaries (div must be zero behind the open boundary) 
    110          !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column 
    111          IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1)   = 0.e0    ! east 
    112          IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1)   = 0.e0    ! west 
    113          IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0    ! north 
    114          IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0    ! south 
    115 #endif 
    116  
    117 #if defined key_bdy 
    118          jgrd=1 !: tracer grid. 
    119          DO jb = 1, nblenrim(jgrd) 
    120            ji = nbi(jb,jgrd) 
    121            jj = nbj(jb,jgrd) 
    122            zhdiv(ji,jj) = 0.e0 
     184     ELSE                                         ! Leap-frog time stepping 
     185         DO jj = 1, jpj 
     186            DO ji = 1, jpi 
     187               ! before <-- now filtered 
     188               sshb  (ji,jj) = sshn(ji,jj)   + atfp * ( sshb  (ji,jj) - 2 * sshn  (ji,jj) + ssha  (ji,jj) )    !& 
     189               ! now <-- after 
     190               sshn  (ji,jj) = ssha  (ji,jj) 
     191            END DO 
    123192         END DO 
    124 #endif 
    125  
    126          CALL lbc_lnk( zhdiv, 'T', 1. ) 
    127  
    128          ! Sea surface elevation time stepping 
    129          ! ----------------------------------- 
    130          zssha(:,:) = sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) ) * tmask(:,:,1) 
    131  
    132          ! Vertical velocity computed from bottom 
    133          ! -------------------------------------- 
    134          DO jk = jpkm1, 1, -1 
    135             wn(:,:,jk) = wn(:,:,jk+1) - fse3t(:,:,jk) * hdivn(:,:,jk) & 
    136               &        - ( zssha(:,:) - sshb(:,:) ) * fse3t_b(:,:,jk) * mut(:,:,jk) / z2dt 
    137          END DO 
    138  
    139       ELSE                             ! Fixed volume  
    140  
    141          ! Vertical velocity computed from bottom 
    142          ! -------------------------------------- 
    143          DO jk = jpkm1, 1, -1 
    144             wn(:,:,jk) = wn(:,:,jk+1) - fse3t(:,:,jk) * hdivn(:,:,jk) 
    145          END DO 
    146        
    147       ENDIF  
    148  
    149       IF(ln_ctl)   CALL prt_ctl(tab3d_1=wn, clinfo1=' w**2 -   : ', mask1=wn) 
    150  
    151    END SUBROUTINE wzv 
     193      ENDIF 
     194 
     195      IF(ln_ctl)   CALL prt_ctl(tab2d_1=sshb    , clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 ) 
     196      ! 
     197   END SUBROUTINE ssh_nxt 
    152198 
    153199   !!====================================================================== 

  • branches/dev_004_VVL/NEMO/OPA_SRC/istate.F90

    r1200 r1368  
    156156         CALL bn2( tb, sb, rn2 )    ! before Brunt Vaissala frequency 
    157157 
    158          IF( .NOT. ln_rstart ) CALL wzv( nit000 )  
    159  
    160158      ENDIF 
    161159 
    162       !                                       ! Vertical velocity 
    163       !                                       ! ----------------- 
    164  
    165       IF( .NOT. lk_vvl )    CALL wzv( nit000 )                         ! from horizontal divergence 
    166       ! 
    167160   END SUBROUTINE istate_init 
    168161 

  • branches/dev_004_VVL/NEMO/OPA_SRC/step.F90

    r1359 r1368  
    201201      ENDIF 
    202202 
     203      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     204      ! Computation of diagnostic variables 
     205      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     206      ! N.B. ua, va, ta, sa arrays are used as workspace in this section 
     207      !----------------------------------------------------------------------- 
     208                       CALL div_cur( kstp )                 ! Horizontal divergence & Relative vorticity 
     209      IF( n_cla == 1 ) CALL div_cla( kstp )                 ! Cross Land Advection (Update Hor. divergence) 
     210                       CALL ssh_wzv( kstp       )           ! after ssh & vertical velocity 
     211      IF( lk_vvl )     CALL dom_vvl                         ! vertical mesh at next time step 
     212 
    203213 
    204214      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    … …  
    360370                               CALL dyn_spg( kstp, indic )    ! surface pressure gradient 
    361371                               CALL dyn_nxt( kstp )           ! lateral velocity at next time step 
    362       IF( lk_vvl )             CALL dom_vvl                   ! vertical mesh at next time step 
    363  
    364  
    365       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    366       ! Computation of diagnostic variables 
    367       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    368       ! N.B. ua, va, ta, sa arrays are used as workspace in this section 
    369       !----------------------------------------------------------------------- 
    370                        CALL div_cur( kstp )                 ! Horizontal divergence & Relative vorticity 
    371       IF( n_cla == 1 ) CALL div_cla( kstp )                 ! Cross Land Advection (Update Hor. divergence) 
    372                        CALL wzv( kstp )                     ! Vertical velocity 
     372 
     373                               CALL ssh_nxt( kstp )           ! sea surface height at next time step 
    373374 
    374375      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
Note: See TracChangeset for help on using the changeset viewer.