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 6225 for branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 – NEMO

Ignore:
Timestamp:
2016-01-08T10:35:19+01:00 (8 years ago)
Author:
jamesharle
Message:

Update MPP_BDY_UPDATE branch to be consistent with head of trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

    r4486 r6225  
    1212 
    1313   !!---------------------------------------------------------------------- 
    14    !!   ssh_nxt        : after ssh 
    15    !!   ssh_swp        : filter ans swap the ssh arrays 
    16    !!   wzv            : compute now vertical velocity 
    17    !!---------------------------------------------------------------------- 
    18    USE oce             ! ocean dynamics and tracers variables 
    19    USE dom_oce         ! ocean space and time domain variables  
    20    USE sbc_oce         ! surface boundary condition: ocean 
    21    USE domvvl          ! Variable volume 
    22    USE divcur          ! hor. divergence and curl      (div & cur routines) 
    23    USE iom             ! I/O library 
    24    USE restart         ! only for lrst_oce 
    25    USE in_out_manager  ! I/O manager 
    26    USE prtctl          ! Print control 
    27    USE phycst 
    28    USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
    29    USE lib_mpp         ! MPP library 
    30    USE bdy_oce 
    31    USE bdy_par          
    32    USE bdydyn2d        ! bdy_ssh routine 
    33    USE diaar5, ONLY:   lk_diaar5 
    34    USE iom 
     14   !!   ssh_nxt       : after ssh 
     15   !!   ssh_swp       : filter ans swap the ssh arrays 
     16   !!   wzv           : compute now vertical velocity 
     17   !!---------------------------------------------------------------------- 
     18   USE oce            ! ocean dynamics and tracers variables 
     19   USE dom_oce        ! ocean space and time domain variables  
     20   USE sbc_oce        ! surface boundary condition: ocean 
     21   USE domvvl         ! Variable volume 
     22   USE divhor         ! horizontal divergence 
     23   USE phycst         ! physical constants 
     24   USE bdy_oce        !  
     25   USE bdy_par        ! 
     26   USE bdydyn2d       ! bdy_ssh routine 
    3527#if defined key_agrif 
    36    USE agrif_opa_update 
    3728   USE agrif_opa_interp 
    3829#endif 
    3930#if defined key_asminc    
    40    USE asminc          ! Assimilation increment 
    41 #endif 
    42    USE wrk_nemo        ! Memory Allocation 
    43    USE timing          ! Timing 
     31   USE   asminc       ! Assimilation increment 
     32#endif 
     33   ! 
     34   USE in_out_manager ! I/O manager 
     35   USE restart        ! only for lrst_oce 
     36   USE prtctl         ! Print control 
     37   USE lbclnk         ! ocean lateral boundary condition (or mpp link) 
     38   USE lib_mpp        ! MPP library 
     39   USE wrk_nemo       ! Memory Allocation 
     40   USE timing         ! Timing 
     41   USE wet_dry         ! Wetting/Drying flux limting 
    4442 
    4543   IMPLICIT NONE 
     
    5149 
    5250   !! * Substitutions 
    53 #  include "domzgr_substitute.h90" 
    5451#  include "vectopt_loop_substitute.h90" 
    5552   !!---------------------------------------------------------------------- 
     
    7067      !!      by the time step. 
    7168      !! 
    72       !! ** action  :   ssha    : after sea surface height 
     69      !! ** action  :   ssha, after sea surface height 
    7370      !! 
    7471      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
    7572      !!---------------------------------------------------------------------- 
    76       ! 
    77       REAL(wp), POINTER, DIMENSION(:,:  ) ::  zhdiv 
    78       INTEGER, INTENT(in) ::   kt                      ! time step 
     73      INTEGER, INTENT(in) ::   kt   ! time step 
    7974      !  
    80       INTEGER             ::   jk                      ! dummy loop indice 
    81       REAL(wp)            ::   z2dt, z1_rau0           ! local scalars 
    82       !!---------------------------------------------------------------------- 
    83       ! 
    84       IF( nn_timing == 1 )  CALL timing_start('ssh_nxt') 
    85       ! 
    86       CALL wrk_alloc( jpi, jpj, zhdiv )  
     75      INTEGER  ::   jk            ! dummy loop indice 
     76      REAL(wp) ::   z2dt, zcoef   ! local scalars 
     77      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhdiv   ! 2D workspace 
     78      !!---------------------------------------------------------------------- 
     79      ! 
     80      IF( nn_timing == 1 )   CALL timing_start('ssh_nxt') 
     81      ! 
     82      CALL wrk_alloc( jpi,jpj,   zhdiv )  
    8783      ! 
    8884      IF( kt == nit000 ) THEN 
    89          ! 
    9085         IF(lwp) WRITE(numout,*) 
    9186         IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height' 
    9287         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
    93          ! 
    94       ENDIF 
    95       ! 
    96       CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity 
     88      ENDIF 
     89      ! 
     90      CALL div_hor( kt )                              ! Horizontal divergence 
    9791      ! 
    9892      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog) 
     
    10498      zhdiv(:,:) = 0._wp 
    10599      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports 
    106         zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk) 
     100        zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk) 
    107101      END DO 
    108102      !                                                ! Sea surface elevation time stepping 
     
    110104      ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 
    111105      !  
    112       z1_rau0 = 0.5_wp * r1_rau0 
    113       ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1) 
    114  
    115 #if ! defined key_dynspg_ts 
    116       ! These lines are not necessary with time splitting since 
    117       ! boundary condition on sea level is set during ts loop 
    118 #if defined key_agrif 
    119       CALL agrif_ssh( kt ) 
    120 #endif 
    121 #if defined key_bdy 
    122       IF (lk_bdy) THEN 
    123          CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary 
    124          CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries 
    125       ENDIF 
    126 #endif 
    127 #endif 
     106      zcoef = 0.5_wp * r1_rau0 
     107 
     108      IF(ln_wd) CALL wad_lmt(sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt) 
     109 
     110      ssha(:,:) = (  sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:) 
     111 
     112      IF ( .NOT.ln_dynspg_ts ) THEN 
     113         ! These lines are not necessary with time splitting since 
     114         ! boundary condition on sea level is set during ts loop 
     115# if defined key_agrif 
     116         CALL agrif_ssh( kt ) 
     117# endif 
     118# if defined key_bdy 
     119         IF( lk_bdy ) THEN 
     120            CALL lbc_lnk( ssha, 'T', 1. )    ! Not sure that's necessary 
     121            CALL bdy_ssh( ssha )             ! Duplicate sea level across open boundaries 
     122         ENDIF 
     123# endif 
     124      ENDIF 
    128125 
    129126#if defined key_asminc 
    130       !                                                ! Include the IAU weighted SSH increment 
    131       IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 
     127      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN     ! Include the IAU weighted SSH increment 
    132128         CALL ssh_asm_inc( kt ) 
    133129         ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:) 
    134130      ENDIF 
    135131#endif 
    136  
    137132      !                                           !------------------------------! 
    138133      !                                           !           outputs            ! 
    139134      !                                           !------------------------------! 
    140       CALL iom_put( "ssh" , sshn                  )   ! sea surface height 
    141       CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height 
    142135      ! 
    143136      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 ) 
     
    165158      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
    166159      !!---------------------------------------------------------------------- 
    167       ! 
    168       INTEGER, INTENT(in) ::   kt           ! time step 
     160      INTEGER, INTENT(in) ::   kt   ! time step 
     161      ! 
     162      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     163      REAL(wp) ::   z1_2dt       ! local scalars 
    169164      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d 
    170165      REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv 
    171       ! 
    172       INTEGER             ::   ji, jj, jk   ! dummy loop indices 
    173       REAL(wp)            ::   z1_2dt       ! local scalars 
    174       !!---------------------------------------------------------------------- 
    175        
    176       IF( nn_timing == 1 )  CALL timing_start('wzv') 
     166      !!---------------------------------------------------------------------- 
     167      ! 
     168      IF( nn_timing == 1 )   CALL timing_start('wzv') 
    177169      ! 
    178170      IF( kt == nit000 ) THEN 
    179          ! 
    180171         IF(lwp) WRITE(numout,*) 
    181172         IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity ' 
     
    183174         ! 
    184175         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all) 
    185          ! 
    186176      ENDIF 
    187177      !                                           !------------------------------! 
     
    199189            DO jj = 2, jpjm1 
    200190               DO ji = fs_2, fs_jpim1   ! vector opt. 
    201                   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) ) 
     191                  zhdiv(ji,jj,jk) = r1_e1e2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) ) 
    202192               END DO 
    203193            END DO 
     
    208198         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence 
    209199            ! computation of w 
    210             wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)                    & 
    211                &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk) 
     200            wn(:,:,jk) = wn(:,:,jk+1) - (  e3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)    & 
     201               &                         + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) )    ) * tmask(:,:,jk) 
    212202         END DO 
    213203         !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0 
     
    216206         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence 
    217207            ! computation of w 
    218             wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk)                                   & 
    219                &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk) 
     208            wn(:,:,jk) = wn(:,:,jk+1) - (  e3t_n(:,:,jk) * hdivn(:,:,jk)                 & 
     209               &                         + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) ) ) * tmask(:,:,jk) 
    220210         END DO 
    221211      ENDIF 
    222212 
    223213#if defined key_bdy 
    224       IF (lk_bdy) THEN 
     214      IF( lk_bdy ) THEN 
    225215         DO jk = 1, jpkm1 
    226216            wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:) 
     
    229219#endif 
    230220      ! 
    231       !                                           !------------------------------! 
    232       !                                           !           outputs            ! 
    233       !                                           !------------------------------! 
    234       CALL iom_put( "woce", wn )   ! vertical velocity 
    235       IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value 
    236          CALL wrk_alloc( jpi, jpj, z2d )  
    237          CALL wrk_alloc( jpi, jpj, jpk, z3d )  
    238          ! Caution: in the VVL case, it only correponds to the baroclinic mass transport. 
    239          z2d(:,:) = rau0 * e12t(:,:) 
    240          DO jk = 1, jpk 
    241             z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:) 
    242          END DO 
    243          CALL iom_put( "w_masstr" , z3d                     )   
    244          CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) 
    245          CALL wrk_dealloc( jpi, jpj, z2d  )  
    246          CALL wrk_dealloc( jpi, jpj, jpk, z3d )  
    247       ENDIF 
    248       ! 
    249221      IF( nn_timing == 1 )  CALL timing_stop('wzv') 
    250  
    251  
     222      ! 
    252223   END SUBROUTINE wzv 
     224 
    253225 
    254226   SUBROUTINE ssh_swp( kt ) 
     
    272244      !!---------------------------------------------------------------------- 
    273245      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     246      ! 
     247      REAL(wp) ::   zcoef   ! local scalar 
    274248      !!---------------------------------------------------------------------- 
    275249      ! 
     
    281255         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
    282256      ENDIF 
    283  
    284 # if defined key_dynspg_ts 
    285       IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ln_bt_fw ) THEN    !** Euler time-stepping: no filter 
    286 # else 
    287       IF ( neuler == 0 .AND. kt == nit000 ) THEN   !** Euler time-stepping at first time-step : no filter 
    288 #endif 
    289          sshb(:,:) = sshn(:,:)                           ! before <-- now 
    290          sshn(:,:) = ssha(:,:)                           ! now    <-- after  (before already = now) 
    291       ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap 
    292          sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )     ! before <-- now filtered 
    293          IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 
    294          sshn(:,:) = ssha(:,:)                           ! now <-- after 
    295       ENDIF 
    296       ! 
    297       ! Update velocity at AGRIF zoom boundaries 
    298 #if defined key_agrif 
    299       IF ( .NOT.Agrif_Root() ) CALL Agrif_Update_Dyn( kt ) 
    300 #endif 
     257      !              !==  Euler time-stepping: no filter, just swap  ==! 
     258      IF(  ( neuler == 0 .AND. kt == nit000 ) .OR.    & 
     259         & ( ln_bt_fw    .AND. ln_dynspg_ts )      ) THEN  
     260         sshb(:,:) = sshn(:,:)                              ! before <-- now 
     261         sshn(:,:) = ssha(:,:)                              ! now    <-- after  (before already = now) 
     262         ! 
     263      ELSE           !==  Leap-Frog time-stepping: Asselin filter + swap  ==! 
     264         !                                                  ! before <-- now filtered 
     265         sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) ) 
     266         IF( .NOT.ln_linssh ) THEN                          ! before <-- with forcing removed 
     267            zcoef = atfp * rdt * r1_rau0 
     268            sshb(:,:) = sshb(:,:) - zcoef * (     emp_b(:,:) - emp   (:,:)   & 
     269               &                             -    rnf_b(:,:) + rnf   (:,:)   & 
     270               &                             + fwfisf_b(:,:) - fwfisf(:,:)   ) * ssmask(:,:) 
     271         ENDIF 
     272         sshn(:,:) = ssha(:,:)                              ! now <-- after 
     273      ENDIF 
    301274      ! 
    302275      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 ) 
    303276      ! 
    304       IF( nn_timing == 1 )  CALL timing_stop('ssh_swp') 
     277      IF( nn_timing == 1 )   CALL timing_stop('ssh_swp') 
    305278      ! 
    306279   END SUBROUTINE ssh_swp 
Note: See TracChangeset for help on using the changeset viewer.