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 6020 for branches/UKMO/icebergs_restart_single_file/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 – NEMO

Ignore:
Timestamp:
2015-12-08T12:39:53+01:00 (8 years ago)
Author:
timgraham
Message:

Merged with head of trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/icebergs_restart_single_file/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90

    r6019 r6020  
    1919   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes 
    2020   !!            3.7  !  2014-04  (G. Madec) add the diagnostic of the time filter trends 
     21   !!            3.7  !  2015-11  (J. Chanut) Free surface simplification 
    2122   !!------------------------------------------------------------------------- 
    2223   
     
    2829   USE sbc_oce         ! Surface boundary condition: ocean fields 
    2930   USE phycst          ! physical constants 
    30    USE dynspg_oce      ! type of surface pressure gradient 
    3131   USE dynadv          ! dynamics: vector invariant versus flux form 
    3232   USE domvvl          ! variable volume 
     
    6868      !!                  ***  ROUTINE dyn_nxt  *** 
    6969      !!                    
    70       !! ** Purpose :   Compute the after horizontal velocity. Apply the boundary  
    71       !!             condition on the after velocity, achieved the time stepping  
     70      !! ** Purpose :   Finalize after horizontal velocity. Apply the boundary  
     71      !!             condition on the after velocity, achieve the time stepping  
    7272      !!             by applying the Asselin filter on now fields and swapping  
    7373      !!             the fields. 
    7474      !! 
    75       !! ** Method  : * After velocity is compute using a leap-frog scheme: 
    76       !!                       (ua,va) = (ub,vb) + 2 rdt (ua,va) 
    77       !!             Note that with flux form advection and variable volume layer 
    78       !!             (lk_vvl=T), the leap-frog is applied on thickness weighted 
    79       !!             velocity. 
    80       !!             Note also that in filtered free surface (lk_dynspg_flt=T), 
    81       !!             the time stepping has already been done in dynspg module 
     75      !! ** Method  : * Ensure after velocities transport matches time splitting 
     76      !!             estimate (ln_dynspg_ts=T) 
    8277      !! 
    8378      !!              * Apply lateral boundary conditions on after velocity  
     
    10196      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    10297      INTEGER  ::   iku, ikv     ! local integers 
    103 #if ! defined key_dynspg_flt 
    104       REAL(wp) ::   z2dt         ! temporary scalar 
    105 #endif 
    106       REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec      ! local scalars 
     98      REAL(wp) ::   zue3a, zue3n, zue3b, zuf           ! local scalars 
    10799      REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      - 
    108100      REAL(wp), POINTER, DIMENSION(:,:)   ::  zue, zve 
     
    112104      IF( nn_timing == 1 )   CALL timing_start('dyn_nxt') 
    113105      ! 
    114       CALL wrk_alloc( jpi,jpj,jpk,  ze3u_f, ze3v_f, zua, zva ) 
    115       IF( lk_dynspg_ts )   CALL wrk_alloc( jpi,jpj, zue, zve ) 
     106      IF( ln_dynspg_ts )   CALL wrk_alloc( jpi,jpj, zue, zve ) 
     107      IF( lk_vvl.AND.(.NOT.ln_dynadv_vec ) ) CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 
     108      IF( l_trddyn     )   CALL wrk_alloc( jpi,jpj,jpk, zua, zva ) 
    116109      ! 
    117110      IF( kt == nit000 ) THEN 
     
    121114      ENDIF 
    122115 
    123 #if defined key_dynspg_flt 
    124       ! 
    125       ! Next velocity :   Leap-frog time stepping already done in dynspg_flt.F routine 
    126       ! ------------- 
    127  
    128       ! Update after velocity on domain lateral boundaries      (only local domain required) 
    129       ! -------------------------------------------------- 
    130       CALL lbc_lnk( ua, 'U', -1. )         ! local domain boundaries 
    131       CALL lbc_lnk( va, 'V', -1. )  
    132       ! 
    133 #else 
    134  
    135 # if defined key_dynspg_exp 
    136       ! Next velocity :   Leap-frog time stepping 
    137       ! ------------- 
    138       z2dt = 2. * rdt                                 ! Euler or leap-frog time step  
    139       IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt 
    140       ! 
    141       IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity 
     116      IF ( ln_dynspg_ts ) THEN 
     117         ! Ensure below that barotropic velocities match time splitting estimate 
     118         ! Compute actual transport and replace it with ts estimate at "after" time step 
     119         zue(:,:) = fse3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1) 
     120         zve(:,:) = fse3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1) 
     121         DO jk = 2, jpkm1 
     122            zue(:,:) = zue(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 
     123            zve(:,:) = zve(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 
     124         END DO 
    142125         DO jk = 1, jpkm1 
    143             ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk) 
    144             va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk) 
    145          END DO 
    146       ELSE                                            ! applied on thickness weighted velocity 
    147          DO jk = 1, jpkm1 
    148             ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      & 
    149                &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   & 
    150                &         / fse3u_a(:,:,jk) * umask(:,:,jk) 
    151             va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      & 
    152                &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   & 
    153                &         / fse3v_a(:,:,jk) * vmask(:,:,jk) 
    154          END DO 
    155       ENDIF 
    156 # endif 
    157  
    158 # if defined key_dynspg_ts 
    159 !!gm IF ( lk_dynspg_ts ) THEN .... 
    160       ! Ensure below that barotropic velocities match time splitting estimate 
    161       ! Compute actual transport and replace it with ts estimate at "after" time step 
    162       zue(:,:) = fse3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1) 
    163       zve(:,:) = fse3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1) 
    164       DO jk = 2, jpkm1 
    165          zue(:,:) = zue(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 
    166          zve(:,:) = zve(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 
    167       END DO 
    168       DO jk = 1, jpkm1 
    169          ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk) 
    170          va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk) 
    171       END DO 
    172  
    173       IF (lk_dynspg_ts.AND.(.NOT.ln_bt_fw)) THEN 
    174          ! Remove advective velocity from "now velocities"  
    175          ! prior to asselin filtering      
    176          ! In the forward case, this is done below after asselin filtering    
    177          ! so that asselin contribution is removed at the same time  
    178          DO jk = 1, jpkm1 
    179             un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk) 
    180             vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk) 
    181          END DO   
    182       ENDIF 
    183 !!gm ENDIF 
    184 # endif 
     126            ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk) 
     127            va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk) 
     128         END DO 
     129 
     130         IF ( .NOT.ln_bt_fw ) THEN 
     131            ! Remove advective velocity from "now velocities"  
     132            ! prior to asselin filtering      
     133            ! In the forward case, this is done below after asselin filtering    
     134            ! so that asselin contribution is removed at the same time  
     135            DO jk = 1, jpkm1 
     136               un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk) 
     137               vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk) 
     138            END DO   
     139         ENDIF 
     140      ENDIF 
    185141 
    186142      ! Update after velocity on domain lateral boundaries 
    187143      ! --------------------------------------------------       
    188       CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries 
    189       CALL lbc_lnk( va, 'V', -1. )  
    190       ! 
    191 # if defined key_bdy 
    192       !                                !* BDY open boundaries 
    193       IF( lk_bdy .AND. lk_dynspg_exp ) CALL bdy_dyn( kt ) 
    194       IF( lk_bdy .AND. lk_dynspg_ts  ) CALL bdy_dyn( kt, dyn3d_only=.true. ) 
    195  
    196 !!$   Do we need a call to bdy_vol here?? 
    197       ! 
    198 # endif 
    199       ! 
    200144# if defined key_agrif 
    201145      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries 
    202146# endif 
    203 #endif 
    204  
     147      ! 
     148      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries 
     149      CALL lbc_lnk( va, 'V', -1. )  
     150      ! 
     151# if defined key_bdy 
     152      !                                !* BDY open boundaries 
     153      IF( lk_bdy .AND. ln_dynspg_exp ) CALL bdy_dyn( kt ) 
     154      IF( lk_bdy .AND. ln_dynspg_ts  ) CALL bdy_dyn( kt, dyn3d_only=.true. ) 
     155 
     156!!$   Do we need a call to bdy_vol here?? 
     157      ! 
     158# endif 
     159      ! 
    205160      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics 
    206161         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step  
     
    259214            ! (used as a now filtered scale factor until the swap) 
    260215            ! ---------------------------------------------------- 
    261             IF (lk_dynspg_ts.AND.ln_bt_fw) THEN 
     216            IF (ln_dynspg_ts.AND.ln_bt_fw) THEN 
    262217               ! No asselin filtering on thicknesses if forward time splitting 
    263                   fse3t_b(:,:,:) = fse3t_n(:,:,:) 
     218               fse3t_b(:,:,:) = fse3t_n(:,:,:) 
    264219            ELSE 
    265220               fse3t_b(:,:,:) = fse3t_n(:,:,:) + atfp * ( fse3t_b(:,:,:) - 2._wp * fse3t_n(:,:,:) + fse3t_a(:,:,:) ) 
    266221               ! Add volume filter correction: compatibility with tracer advection scheme 
    267222               ! => time filter + conservation correction (only at the first level) 
    268                fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) & 
    269                               &                                          -rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1) 
     223               IF ( nn_isf == 0) THEN   ! if no ice shelf melting 
     224                  fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) & 
     225                                 &                                          -rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1) 
     226               ELSE                     ! if ice shelf melting 
     227                  DO jj = 1,jpj 
     228                     DO ji = 1,jpi 
     229                        jk = mikt(ji,jj) 
     230                        fse3t_b(ji,jj,jk) = fse3t_b(ji,jj,jk) - atfp * rdt * r1_rau0                       & 
     231                                          &                          * ( (emp_b(ji,jj)    - emp(ji,jj)   ) & 
     232                                          &                            - (rnf_b(ji,jj)    - rnf(ji,jj)   ) & 
     233                                          &                            + (fwfisf_b(ji,jj) - fwfisf(ji,jj)) ) * tmask(ji,jj,jk) 
     234                     END DO 
     235                  END DO 
     236               END IF 
    270237            ENDIF 
    271238            ! 
     
    324291         ENDIF 
    325292         ! 
    326          IF (lk_dynspg_ts.AND.ln_bt_fw) THEN 
     293         IF (ln_dynspg_ts.AND.ln_bt_fw) THEN 
    327294            ! Revert "before" velocities to time split estimate 
    328295            ! Doing it here also means that asselin filter contribution is removed   
     
    388355      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   & 
    389356         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask ) 
    390       !  
    391       CALL wrk_dealloc( jpi,jpj,jpk,  ze3u_f, ze3v_f, zua, zva ) 
    392       IF( lk_dynspg_ts )   CALL wrk_dealloc( jpi,jpj, zue, zve ) 
     357      ! 
     358      IF( ln_dynspg_ts )   CALL wrk_dealloc( jpi,jpj, zue, zve ) 
     359      IF( lk_vvl.AND.(.NOT.ln_dynadv_vec ) ) CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 
     360      IF( l_trddyn     )   CALL wrk_dealloc( jpi,jpj,jpk, zua, zva ) 
    393361      ! 
    394362      IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt') 
Note: See TracChangeset for help on using the changeset viewer.