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 7422 for branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Ignore:
Timestamp:
2016-12-01T18:17:41+01:00 (8 years ago)
Author:
gm
Message:

#1805 dev_INGV_UKMO_2016: improve Stokes drift (including dynspg_ts , Stokes-Coriolis force , and GLS surface roughness

Location:
branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/DYN
Files:
1 deleted
2 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r7351 r7422  
    1515   !!             3.6  ! 2013-11  (A. Coward) Update for z-tilde compatibility 
    1616   !!             3.7  ! 2015-11  (J. Chanut) free surface simplification 
     17   !!              -   ! 2016-12  (G. Madec, E. Clementi) update for Stoke-Drift divergence 
    1718   !!--------------------------------------------------------------------- 
    1819 
     
    3839   USE sbctide         ! tides 
    3940   USE updtide         ! tide potential 
     41   USE sbcwave         ! surface wave 
    4042   ! 
    4143   USE in_out_manager  ! I/O manager 
     
    168170      CALL wrk_alloc( jpi,jpj,   zsshu_a, zsshv_a                  ) 
    169171      CALL wrk_alloc( jpi,jpj,   zhf ) 
    170       IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 ) 
    171       ! 
    172       zmdi=1.e+20                               !  missing data indicator for masking 
     172      IF( ln_wd )   CALL wrk_alloc( jpi,jpj,   zcpx, zcpy, wduflt1, wdvflt1 ) 
     173      ! 
    173174      !                                         !* Local constant initialization 
    174       z1_12 = 1._wp / 12._wp  
     175      zmdi = 1.e+20                                   !  missing data indicator for masking 
     176      ! 
     177      z1_12 = 1._wp / 12._wp                          ! constants 
    175178      z1_8  = 0.125_wp                                    
    176179      z1_4  = 0.25_wp 
    177180      z1_2  = 0.5_wp      
    178181      zraur = 1._wp / rau0 
    179       !                                            ! reciprocal of baroclinic time step  
     182      !                                               ! baroclinic time step  
    180183      IF( kt == nit000 .AND. neuler == 0 ) THEN   ;   z2dt_bf =          rdt 
    181184      ELSE                                        ;   z2dt_bf = 2.0_wp * rdt 
    182185      ENDIF 
    183       z1_2dt_b = 1.0_wp / z2dt_bf  
    184       ! 
    185       ll_init     = ln_bt_av                       ! if no time averaging, then no specific restart  
     186      z1_2dt_b = 1.0_wp / z2dt_bf                     ! reciprocal of baroclinic time step  
     187      ! 
     188      ll_init     = ln_bt_av                          ! if no time averaging, then no specific restart  
    186189      ll_fw_start = .FALSE. 
    187       !                                            ! time offset in steps for bdy data update 
     190      !                                               ! time offset in steps for bdy data update 
    188191      IF( .NOT.ln_bt_fw ) THEN   ;   noffset = - nn_baro 
    189192      ELSE                       ;   noffset =   0  
     
    255258            zwz(:,:) = 0._wp 
    256259            zhf(:,:) = 0._wp 
    257             IF ( .not. ln_sco ) THEN 
     260            IF( .not. ln_sco ) THEN 
    258261 
    259262!!gm  agree the JC comment  : this should be done in a much clear way 
     
    324327         END DO 
    325328      END DO 
     329       
     330!!gm  Question here when removing the Vertically integrated trends, we remove the vertically integrated NL trends on momentum.... 
     331!!gm  Is it correct to do so ?   I think so... 
     332       
     333       
    326334      !                                   !* barotropic Coriolis trends (vorticity scheme dependent) 
    327335      !                                   ! -------------------------------------------------------- 
     
    373381      !                                   ! ---------------------------------------------------- 
    374382      IF( .NOT.ln_linssh ) THEN                 ! Variable volume : remove surface pressure gradient 
    375         IF( ln_wd ) THEN                        ! Calculating and applying W/D gravity filters 
    376           wduflt1(:,:) = 1.0_wp 
    377           wdvflt1(:,:) = 1.0_wp 
    378           DO jj = 2, jpjm1 
    379              DO ji = 2, jpim1 
    380                 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))   & 
     383         IF( ln_wd ) THEN                        ! Calculating and applying W/D gravity filters 
     384            wduflt1(:,:) = 1.0_wp 
     385            wdvflt1(:,:) = 1.0_wp 
     386            DO jj = 2, jpjm1 
     387               DO ji = 2, jpim1 
     388                  ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))   & 
    381389                        & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj))   & 
    382390                        &  > rn_wdmin1 + rn_wdmin2 
    383                 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))   & 
     391                  ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))   & 
    384392                        &                                   + rn_wdmin1 + rn_wdmin2 
    385                 IF(ll_tmp1) THEN 
    386                   zcpx(ji,jj)    = 1.0_wp 
    387                 ELSEIF(ll_tmp2) THEN 
    388                    ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen here 
    389                   zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) & 
     393                  IF(ll_tmp1) THEN 
     394                     zcpx(ji,jj)    = 1.0_wp 
     395                  ELSEIF(ll_tmp2) THEN 
     396                     ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen here 
     397                     zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) & 
    390398                        &          /(sshn(ji+1,jj) - sshn(ji,jj))) 
    391                 ELSE 
    392                   zcpx(ji,jj)    = 0._wp 
    393                   wduflt1(ji,jj) = 0.0_wp 
    394                 END IF 
    395  
    396                 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))   & 
     399                  ELSE 
     400                     zcpx(ji,jj)    = 0._wp 
     401                     wduflt1(ji,jj) = 0.0_wp 
     402                  END IF 
     403 
     404                  ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))   & 
    397405                        & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1))   & 
    398406                        &  > rn_wdmin1 + rn_wdmin2 
    399                 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))   & 
     407                  ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))   & 
    400408                        &                                   + rn_wdmin1 + rn_wdmin2 
    401                 IF(ll_tmp1) THEN 
    402                    zcpy(ji,jj)    = 1.0_wp 
    403                 ELSEIF(ll_tmp2) THEN 
    404                    ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen here 
    405                   zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) & 
     409                  IF(ll_tmp1) THEN 
     410                     zcpy(ji,jj)    = 1.0_wp 
     411                  ELSEIF(ll_tmp2) THEN 
     412                     ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen here 
     413                     zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) & 
    406414                        &          /(sshn(ji,jj+1) - sshn(ji,jj))) 
    407                 ELSE 
    408                   zcpy(ji,jj)    = 0._wp 
    409                   wdvflt1(ji,jj) = 0.0_wp 
    410                 ENDIF 
    411  
    412              END DO 
    413            END DO 
    414  
    415            CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
    416  
    417            DO jj = 2, jpjm1 
    418               DO ji = 2, jpim1 
    419                  zu_trd(ji,jj) = ( zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   & 
    420                         &                        * r1_e1u(ji,jj) ) * zcpx(ji,jj) * wduflt1(ji,jj) 
    421                  zv_trd(ji,jj) = ( zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) )   & 
    422                         &                        * r1_e2v(ji,jj) ) * zcpy(ji,jj) * wdvflt1(ji,jj) 
    423               END DO 
    424            END DO 
    425  
    426          ELSE 
    427  
    428            DO jj = 2, jpjm1 
    429               DO ji = fs_2, fs_jpim1   ! vector opt. 
    430                  zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) * r1_e1u(ji,jj) 
    431                  zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) * r1_e2v(ji,jj)  
    432               END DO 
    433            END DO 
    434         ENDIF 
    435  
    436       ENDIF 
     415                  ELSE 
     416                     zcpy(ji,jj)    = 0._wp 
     417                     wdvflt1(ji,jj) = 0.0_wp 
     418                  ENDIF 
     419 
     420               END DO 
     421            END DO 
     422            ! 
     423            CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     424            ! 
     425            DO jj = 2, jpjm1 
     426               DO ji = 2, jpim1 
     427                  zu_trd(ji,jj) = ( zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   & 
     428                         &                        * r1_e1u(ji,jj) ) * zcpx(ji,jj) * wduflt1(ji,jj) 
     429                  zv_trd(ji,jj) = ( zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) )   & 
     430                         &                        * r1_e2v(ji,jj) ) * zcpy(ji,jj) * wdvflt1(ji,jj) 
     431               END DO 
     432            END DO 
     433            ! 
     434         ELSE        ! no Wet & Drying 
     435            ! 
     436            DO jj = 2, jpjm1 
     437               DO ji = fs_2, fs_jpim1   ! vector opt. 
     438                  zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) * r1_e1u(ji,jj) 
     439                  zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) * r1_e2v(ji,jj)  
     440               END DO 
     441            END DO 
     442         ENDIF 
     443         ! 
     444      ENDIF          ! end non linear free surface case 
    437445 
    438446      DO jj = 2, jpjm1                          ! Remove coriolis term (and possibly spg) from barotropic trend 
     
    473481      END IF 
    474482      ! 
     483!!gm   TOP stress only !!!   this should be with a test on ISF use or not 
    475484      !                                         ! Add top stress contribution from baroclinic velocities:       
    476       IF (ln_bt_fw) THEN 
     485      IF( ln_bt_fw ) THEN 
    477486         DO jj = 2, jpjm1 
    478487            DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    538547                &                        + fwfisf(:,:) + fwfisf_b(:,:)                     ) 
    539548      ENDIF 
     549      ! 
     550      IF( ln_sdw ) THEN                         ! Stokes drift divergence added if necessary 
     551         zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:) 
     552      ENDIF 
     553      ! 
    540554#if defined key_asminc 
    541555      !                                         ! Include the IAU weighted SSH increment 
  • branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90

    r7351 r7422  
    1717   !!            3.7  ! 2014-04  (G. Madec) trend simplification: suppress jpdyn_trd_dat vorticity  
    1818   !!             -   ! 2014-06  (G. Madec) suppression of velocity curl from in-core memory 
     19   !!             -   ! 2016-12  (G. Madec, E. Clementi) add Stokes-Coriolis trends (ln_stcor=T) 
    1920   !!---------------------------------------------------------------------- 
    2021 
     
    3233   USE trd_oce        ! trends: ocean variables 
    3334   USE trddyn         ! trend manager: dynamics 
     35   USE sbcwave        ! Surface Waves (add Stokes-Coriolis force) 
     36   USE sbc_oce , ONLY : ln_stcor    ! use Stoke-Coriolis force 
    3437   ! 
    3538   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     
    7780#  include "vectopt_loop_substitute.h90" 
    7881   !!---------------------------------------------------------------------- 
    79    !! NEMO/OPA 3.7 , NEMO Consortium (2014) 
     82   !! NEMO/OPA 3.7 , NEMO Consortium (2016) 
    8083   !! $Id$ 
    8184   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    108111            ztrdu(:,:,:) = ua(:,:,:) 
    109112            ztrdv(:,:,:) = va(:,:,:) 
    110             CALL vor_ene( kt, nrvm, ua, va )                      ! relative vorticity or metric trend 
     113            CALL vor_ene( kt, nrvm, un , vn , ua, va )                    ! relative vorticity or metric trend 
    111114            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    112115            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     
    114117            ztrdu(:,:,:) = ua(:,:,:) 
    115118            ztrdv(:,:,:) = va(:,:,:) 
    116             CALL vor_ene( kt, ncor, ua, va )                      ! planetary vorticity trend 
     119            CALL vor_ene( kt, ncor, un , vn , ua, va )                    ! planetary vorticity trend 
    117120            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    118121            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    119122            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 
    120          ELSE 
    121             CALL vor_ene( kt, ntot, ua, va )                ! total vorticity trend 
     123         ELSE                                               ! total vorticity trend 
     124                             CALL vor_ene( kt, ntot, un , vn , ua, va )   ! total vorticity trend 
     125            IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend 
    122126         ENDIF 
    123127         ! 
     
    126130            ztrdu(:,:,:) = ua(:,:,:) 
    127131            ztrdv(:,:,:) = va(:,:,:) 
    128             CALL vor_ens( kt, nrvm, ua, va )                      ! relative vorticity or metric trend 
     132            CALL vor_ens( kt, nrvm, un , vn , ua, va )            ! relative vorticity or metric trend 
    129133            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    130134            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     
    132136            ztrdu(:,:,:) = ua(:,:,:) 
    133137            ztrdv(:,:,:) = va(:,:,:) 
    134             CALL vor_ens( kt, ncor, ua, va )                      ! planetary vorticity trend 
     138            CALL vor_ens( kt, ncor, un , vn , ua, va )            ! planetary vorticity trend 
    135139            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    136140            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    137141            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 
    138          ELSE 
    139             CALL vor_ens( kt, ntot, ua, va )                ! total vorticity trend 
     142         ELSE                                               ! total vorticity trend 
     143                             CALL vor_ens( kt, ntot, un , vn , ua, va )  ! total vorticity trend 
     144            IF( ln_stcor )   CALL vor_ens( kt, ncor, usd, vsd, ua, va )  ! add the Stokes-Coriolis trend 
    140145         ENDIF 
    141146         ! 
     
    144149            ztrdu(:,:,:) = ua(:,:,:) 
    145150            ztrdv(:,:,:) = va(:,:,:) 
    146             CALL vor_ens( kt, nrvm, ua, va )                      ! relative vorticity or metric trend (ens) 
     151            CALL vor_ens( kt, nrvm, un , vn , ua, va )            ! relative vorticity or metric trend (ens) 
    147152            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    148153            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     
    150155            ztrdu(:,:,:) = ua(:,:,:) 
    151156            ztrdv(:,:,:) = va(:,:,:) 
    152             CALL vor_ene( kt, ncor, ua, va )                      ! planetary vorticity trend (ene) 
     157            CALL vor_ene( kt, ncor, un , vn , ua, va )            ! planetary vorticity trend (ene) 
    153158            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    154159            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    155160            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 
    156          ELSE 
    157             CALL vor_ens( kt, nrvm, ua, va )                ! relative vorticity or metric trend (ens) 
    158             CALL vor_ene( kt, ncor, ua, va )                ! planetary vorticity trend (ene) 
     161         ELSE                                               ! total vorticity trend 
     162                             CALL vor_ens( kt, nrvm, un , vn , ua, va )   ! relative vorticity or metric trend (ens) 
     163                             CALL vor_ene( kt, ncor, un , vn , ua, va )   ! planetary vorticity trend (ene) 
     164            IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend 
    159165        ENDIF 
    160166         ! 
     
    163169            ztrdu(:,:,:) = ua(:,:,:) 
    164170            ztrdv(:,:,:) = va(:,:,:) 
    165             CALL vor_een( kt, nrvm, ua, va )                      ! relative vorticity or metric trend 
     171            CALL vor_een( kt, nrvm, un , vn , ua, va )            ! relative vorticity or metric trend 
    166172            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    167173            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     
    169175            ztrdu(:,:,:) = ua(:,:,:) 
    170176            ztrdv(:,:,:) = va(:,:,:) 
    171             CALL vor_een( kt, ncor, ua, va )                      ! planetary vorticity trend 
     177            CALL vor_een( kt, ncor, un , vn , ua, va )            ! planetary vorticity trend 
    172178            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    173179            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    174180            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 
    175          ELSE 
    176             CALL vor_een( kt, ntot, ua, va )                ! total vorticity trend 
     181         ELSE                                               ! total vorticity trend 
     182                             CALL vor_een( kt, ntot, un , vn , ua, va )   ! total vorticity trend 
     183            IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend 
    177184         ENDIF 
    178185         ! 
     
    190197 
    191198 
    192    SUBROUTINE vor_ene( kt, kvor, pua, pva ) 
     199   SUBROUTINE vor_ene( kt, kvor, pun, pvn, pua, pva ) 
    193200      !!---------------------------------------------------------------------- 
    194201      !!                  ***  ROUTINE vor_ene  *** 
     
    210217      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 
    211218      !!---------------------------------------------------------------------- 
    212       INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index 
    213       INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ; 
    214       !                                                           ! =nrvm (relative vorticity or metric) 
    215       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend 
    216       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend 
     219      INTEGER , INTENT(in   )                         ::   kt          ! ocean time-step index 
     220      INTEGER , INTENT(in   )                         ::   kvor        ! =ncor (planetary) ; =ntot (total) ; 
     221      !                                                                ! =nrvm (relative vorticity or metric) 
     222      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities 
     223      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua, pva    ! total v-trend 
    217224      ! 
    218225      INTEGER  ::   ji, jj, jk           ! dummy loop indices 
     
    223230      IF( nn_timing == 1 )  CALL timing_start('vor_ene') 
    224231      ! 
    225       CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz )  
     232      CALL wrk_alloc( jpi,jpj,  zwx, zwy, zwz )  
    226233      ! 
    227234      IF( kt == nit000 ) THEN 
     
    241248            DO jj = 1, jpjm1 
    242249               DO ji = 1, fs_jpim1   ! vector opt. 
    243                   zwz(ji,jj) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    & 
    244                      &          - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
     250                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
     251                     &          - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    245252               END DO 
    246253            END DO 
     
    248255            DO jj = 1, jpjm1 
    249256               DO ji = 1, fs_jpim1   ! vector opt. 
    250                   zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    251                        &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
     257                  zwz(ji,jj) = (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
     258                       &         - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    252259                       &     * 0.5 * r1_e1e2f(ji,jj) 
    253260               END DO 
     
    256263            DO jj = 1, jpjm1 
    257264               DO ji = 1, fs_jpim1   ! vector opt. 
    258                   zwz(ji,jj) = ff(ji,jj) + (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    & 
    259                      &                      - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) & 
     265                  zwz(ji,jj) = ff(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
     266                     &                      - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) & 
    260267                     &                   * r1_e1e2f(ji,jj) 
    261268               END DO 
     
    265272               DO ji = 1, fs_jpim1   ! vector opt. 
    266273                  zwz(ji,jj) = ff(ji,jj)                                                                        & 
    267                        &     + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    268                        &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
     274                       &     + (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
     275                       &         - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    269276                       &     * 0.5 * r1_e1e2f(ji,jj) 
    270277               END DO 
     
    282289         ENDIF 
    283290 
    284          IF( ln_sco ) THEN 
    285             zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk) 
    286             zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 
    287             zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 
    288          ELSE 
    289             zwx(:,:) = e2u(:,:) * un(:,:,jk) 
    290             zwy(:,:) = e1v(:,:) * vn(:,:,jk) 
    291          ENDIF 
     291         zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk) 
     292         zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk) 
     293         zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk) 
     294 
    292295         !                                   !==  compute and add the vorticity term trend  =! 
    293296         DO jj = 2, jpjm1 
     
    304307      END DO                                           !   End of slab 
    305308      !                                                ! =============== 
    306       CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz )  
     309      CALL wrk_dealloc( jpi,jpj,  zwx, zwy, zwz )  
    307310      ! 
    308311      IF( nn_timing == 1 )  CALL timing_stop('vor_ene') 
     
    311314 
    312315 
    313    SUBROUTINE vor_ens( kt, kvor, pua, pva ) 
     316   SUBROUTINE vor_ens( kt, kvor, pun, pvn, pua, pva ) 
    314317      !!---------------------------------------------------------------------- 
    315318      !!                ***  ROUTINE vor_ens  *** 
     
    331334      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 
    332335      !!---------------------------------------------------------------------- 
    333       INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index 
    334       INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ; 
    335          !                                                        ! =nrvm (relative vorticity or metric) 
    336       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend 
    337       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend 
     336      INTEGER , INTENT(in   )                         ::   kt          ! ocean time-step index 
     337      INTEGER , INTENT(in   )                         ::   kvor        ! =ncor (planetary) ; =ntot (total) ; 
     338         !                                                             ! =nrvm (relative vorticity or metric) 
     339      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities 
     340      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua, pva    ! total v-trend 
    338341      ! 
    339342      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     
    344347      IF( nn_timing == 1 )  CALL timing_start('vor_ens') 
    345348      ! 
    346       CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz )  
     349      CALL wrk_alloc( jpi,jpj,  zwx, zwy, zwz )  
    347350      ! 
    348351      IF( kt == nit000 ) THEN 
     
    361364            DO jj = 1, jpjm1 
    362365               DO ji = 1, fs_jpim1   ! vector opt. 
    363                   zwz(ji,jj) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    & 
    364                      &          - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
     366                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
     367                     &          - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
    365368               END DO 
    366369            END DO 
     
    368371            DO jj = 1, jpjm1 
    369372               DO ji = 1, fs_jpim1   ! vector opt. 
    370                   zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    371                        &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
     373                  zwz(ji,jj) = (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
     374                       &         - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    372375                       &     * 0.5 * r1_e1e2f(ji,jj) 
    373376               END DO 
     
    376379            DO jj = 1, jpjm1 
    377380               DO ji = 1, fs_jpim1   ! vector opt. 
    378                   zwz(ji,jj) = ff(ji,jj) + (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    & 
    379                      &                      - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) & 
     381                  zwz(ji,jj) = ff(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
     382                     &                      - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) & 
    380383                     &                   * r1_e1e2f(ji,jj) 
    381384               END DO 
     
    385388               DO ji = 1, fs_jpim1   ! vector opt. 
    386389                  zwz(ji,jj) = ff(ji,jj)                                                                       & 
    387                        &     + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    388                        &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
     390                       &     + (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
     391                       &         - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    389392                       &     * 0.5 * r1_e1e2f(ji,jj) 
    390393               END DO 
     
    402405         ENDIF 
    403406         ! 
    404          IF( ln_sco ) THEN                   !==  horizontal fluxes  ==! 
    405             zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk) 
    406             zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 
    407             zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 
    408          ELSE 
    409             zwx(:,:) = e2u(:,:) * un(:,:,jk) 
    410             zwy(:,:) = e1v(:,:) * vn(:,:,jk) 
    411          ENDIF 
     407         !                                 !==  horizontal fluxes  ==! 
     408         zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk) 
     409         zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk) 
     410         zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk) 
     411         ! 
    412412         !                                   !==  compute and add the vorticity term trend  =! 
    413413         DO jj = 2, jpjm1 
     
    424424      END DO                                           !   End of slab 
    425425      !                                                ! =============== 
    426       CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz )  
     426      CALL wrk_dealloc( jpi,jpj,  zwx, zwy, zwz )  
    427427      ! 
    428428      IF( nn_timing == 1 )  CALL timing_stop('vor_ens') 
     
    431431 
    432432 
    433    SUBROUTINE vor_een( kt, kvor, pua, pva ) 
     433   SUBROUTINE vor_een( kt, kvor, pun, pvn, pua, pva ) 
    434434      !!---------------------------------------------------------------------- 
    435435      !!                ***  ROUTINE vor_een  *** 
     
    448448      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36 
    449449      !!---------------------------------------------------------------------- 
    450       INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index 
    451       INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ; =nrvm (relative or metric) 
    452       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend 
    453       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend 
     450      INTEGER , INTENT(in   )                         ::   kt          ! ocean time-step index 
     451      INTEGER , INTENT(in   )                         ::   kvor        ! =ncor (planetary) ; =ntot (total) ; 
     452         !                                                             ! =nrvm (relative vorticity or metric) 
     453      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities 
     454      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua, pva    ! total v-trend 
    454455      ! 
    455456      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     
    512513            DO jj = 1, jpjm1 
    513514               DO ji = 1, fs_jpim1   ! vector opt. 
    514                   zwz(ji,jj) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    & 
    515                      &          - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) & 
     515                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
     516                     &          - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) & 
    516517                     &       * r1_e1e2f(ji,jj) * z1_e3f(ji,jj) 
    517518               END DO 
     
    520521            DO jj = 1, jpjm1 
    521522               DO ji = 1, fs_jpim1   ! vector opt. 
    522                   zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    523                        &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
     523                  zwz(ji,jj) = (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
     524                       &         - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    524525                       &     * 0.5 * r1_e1e2f(ji,jj) * z1_e3f(ji,jj) 
    525526               END DO 
     
    528529            DO jj = 1, jpjm1 
    529530               DO ji = 1, fs_jpim1   ! vector opt. 
    530                   zwz(ji,jj) = (  ff(ji,jj) + (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    & 
    531                      &                         - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) & 
     531                  zwz(ji,jj) = (  ff(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    & 
     532                     &                         - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) & 
    532533                     &                      * r1_e1e2f(ji,jj)    ) * z1_e3f(ji,jj) 
    533534               END DO 
     
    537538               DO ji = 1, fs_jpim1   ! vector opt. 
    538539                  zwz(ji,jj) = (  ff(ji,jj)                                                                        & 
    539                        &        + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    540                        &            - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
     540                       &        + (   ( pvn(ji+1,jj  ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
     541                       &            - ( pun(ji  ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    541542                       &        * 0.5 * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj) 
    542543               END DO 
     
    557558         ! 
    558559         !                                   !==  horizontal fluxes  ==! 
    559          zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 
    560          zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 
     560         zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk) 
     561         zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk) 
    561562 
    562563         !                                   !==  compute and add the vorticity term trend  =! 
Note: See TracChangeset for help on using the changeset viewer.