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

Changeset 7481


Ignore:
Timestamp:
2016-12-08T19:33:38+01:00 (7 years ago)
Author:
jcastill
Message:

Changes as in branches/2016/dev_INGV_UKMO_2016@7451

Location:
branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM
Files:
1 deleted
12 edited

Legend:

Unmodified
Added
Removed

  • branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/CONFIG/SHARED/namelist_ref

    r7471 r7481  
    10041004   rn_hsro       =  0.02   !  Minimum surface roughness 
    10051005   rn_frac_hs    =   1.3   !  Fraction of wave height as roughness (if nn_z0_met=2) 
    1006    nn_z0_met     =     2   !  Method for surface roughness computation (0/1/2) 
     1006   nn_z0_met     =     2   !  Method for surface roughness computation (0/1/2/3) 
     1007      !                             ! =3 requires ln_wave=T 
    10071008   nn_bc_surf    =     1   !  surface condition (0/1=Dir/Neum) 
    10081009   nn_bc_bot     =     1   !  bottom condition (0/1=Dir/Neum) 
    … …  
    12821283/ 
    12831284!----------------------------------------------------------------------- 
    1284 &namsbc_wave   ! External fields from wave model 
     1285&namsbc_wave   ! External fields from wave model                        (ln_wave=T) 
    12851286!----------------------------------------------------------------------- 
    12861287!              !  file name  ! frequency (hours) ! variable     ! time interp. !  clim   ! 'yearly'/ ! weights  ! rotation ! land/sea mask ! 
    … …  
    12891290   sn_usd      =  'sdw_wave' ,        1          , 'u_sd2d'     ,     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
    12901291   sn_vsd      =  'sdw_wave' ,        1          , 'v_sd2d'     ,     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
    1291    sn_swh      =  'sdw_wave' ,        1          , 'hs'         ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''  
     1292   sn_hsw      =  'sdw_wave' ,        1          , 'hs'         ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''  
    12921293   sn_wmp      =  'sdw_wave' ,        1          , 'wmp'        ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''  
    12931294   sn_wnum     =  'sdw_wave' ,        1          , 'wave_num'   ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''  

  • branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r7470 r7481  
    1111   !!             3.5  ! 2013-07  (J. Chanut) Switch to Forward-backward time stepping 
    1212   !!             3.6  ! 2013-11  (A. Coward) Update for z-tilde compatibility 
     13   !!              -   ! 2016-12  (G. Madec, E. Clementi) update for Stoke-Drift divergence 
    1314   !!--------------------------------------------------------------------- 
    1415#if defined key_dynspg_ts   ||   defined key_esopa 
    … …  
    3132   USE sbctide         ! tides 
    3233   USE updtide         ! tide potential 
     34   USE sbcwave         ! surface wave 
     35   ! 
     36   USE sbcwave         ! surface wave 
    3337   USE lib_mpp         ! distributed memory computing library 
    3438   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    … …  
    459463                &                        + fwfisf(:,:) + fwfisf_b(:,:)                     ) 
    460464      ENDIF 
     465      !  
     466      IF( ln_sdw ) THEN                         ! Stokes drift divergence added if necessary  
     467         zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:)  
     468      ENDIF  
     469      !  
    461470#if defined key_asminc 
    462471      !                                         ! Include the IAU weighted SSH increment 

  • branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90

    r7470 r7481  
    3232   USE trd_oce        ! trends: ocean variables 
    3333   USE trddyn         ! trend manager: dynamics 
     34   USE sbcwave        ! Surface Waves (add Stokes-Coriolis force)  
     35   USE sbc_oce , ONLY : ln_stcor    ! use Stoke-Coriolis force  
     36   !  
    3437   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    3538   USE prtctl         ! Print control 
    … …  
    9194      ! 
    9295      CASE ( -1 )                                      ! esopa: test all possibility with control print 
    93          CALL vor_ene( kt, ntot, ua, va ) 
     96         CALL vor_ene( kt, ntot, un, vn, ua, va ) 
    9497         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor0 - Ua: ', mask1=umask, & 
    9598            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    96          CALL vor_ens( kt, ntot, ua, va ) 
     99         CALL vor_ens( kt, ntot, un, vn, ua, va ) 
    97100         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor1 - Ua: ', mask1=umask, & 
    98101            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    … …  
    100103         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor2 - Ua: ', mask1=umask, & 
    101104            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    102          CALL vor_een( kt, ntot, ua, va ) 
     105         CALL vor_een( kt, ntot, un, vn, ua, va ) 
    103106         CALL prt_ctl( tab3d_1=ua, clinfo1=' vor3 - Ua: ', mask1=umask, & 
    104107            &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    … …  
    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 ) 
    120123         ELSE 
    121             CALL vor_ene( kt, ntot, ua, va )                ! total vorticity 
     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 ) 
    138142         ELSE 
    139             CALL vor_ens( kt, ntot, ua, va )                ! total vorticity 
     143                             CALL vor_ens( kt, ntot, un, vn, ua, va )    ! total vorticity 
     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 ) 
    156161         ELSE 
    157             CALL vor_mix( kt )                               ! total vorticity (mix=ens-ene) 
     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 
    158165         ENDIF 
    159166         ! 
    … …  
    162169            ztrdu(:,:,:) = ua(:,:,:) 
    163170            ztrdv(:,:,:) = va(:,:,:) 
    164             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 
    165172            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    166173            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    … …  
    168175            ztrdu(:,:,:) = ua(:,:,:) 
    169176            ztrdv(:,:,:) = va(:,:,:) 
    170             CALL vor_een( kt, ncor, ua, va )                ! planetary vorticity trend 
     177            CALL vor_een( kt, ncor, un, vn, ua, va )        ! planetary vorticity trend 
    171178            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    172179            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    173180            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 
    174181         ELSE 
    175             CALL vor_een( kt, ntot, ua, va )                ! total vorticity 
     182                             CALL vor_een( kt, ntot, un, vn, ua, va )    ! total vorticity 
     183            IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )  ! add the Stokes-Coriolis trend 
    176184         ENDIF 
    177185         ! 
    … …  
    189197 
    190198 
    191    SUBROUTINE vor_ene( kt, kvor, pua, pva ) 
     199   SUBROUTINE vor_ene( kt, kvor, pun, pvn, pua, pva ) 
    192200      !!---------------------------------------------------------------------- 
    193201      !!                  ***  ROUTINE vor_ene  *** 
    … …  
    219227      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend 
    220228      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend 
     229      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities 
    221230      ! 
    222231      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    … …  
    246255         SELECT CASE( kvor )      ! vorticity considered 
    247256         CASE ( 1 )   ;   zwz(:,:) =                  ff(:,:)      ! planetary vorticity (Coriolis) 
    248          CASE ( 2 )   ;   zwz(:,:) =   rotn(:,:,jk)                ! relative  vorticity 
     257         CASE ( 2 )                                                ! relative  vorticity 
     258            DO jj = 1, jpjm1  
     259               DO ji = 1, fs_jpim1   ! vector opt.  
     260                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &  
     261                     &          - e1u(ji  ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e12f(ji,jj)  
     262               END DO  
     263            END DO  
    249264         CASE ( 3 )                                                ! metric term 
    250265            DO jj = 1, jpjm1 
    251266               DO ji = 1, fs_jpim1   ! vector opt. 
    252                   zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    253                        &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    254                        &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) 
     267                  zwz(ji,jj) = (   ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &  
     268                       &         - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
     269                       &     * 0.5 * r1_e12f(ji,jj) 
    255270               END DO 
    256271            END DO 
    257          CASE ( 4 )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) )    ! total (relative + planetary vorticity) 
     272         CASE ( 4 )                                                ! total (relative + planetary vorticity) 
     273            DO jj = 1, jpjm1  
     274               DO ji = 1, fs_jpim1   ! vector opt.  
     275                  zwz(ji,jj) = ff(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &  
     276                     &                      - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) &  
     277                     &                   * r1_e12f(ji,jj)  
     278               END DO  
     279            END DO 
    258280         CASE ( 5 )                                                ! total (coriolis + metric) 
    259             DO jj = 1, jpjm1 
    260                DO ji = 1, fs_jpim1   ! vector opt. 
    261                   zwz(ji,jj) = ( ff (ji,jj)                                                                       & 
    262                        &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    263                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    264                        &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                               & 
    265                        &       ) 
    266                END DO 
     281            DO jj = 1, jpjm1  
     282               DO ji = 1, fs_jpim1   ! vector opt.  
     283                  zwz(ji,jj) = ff(ji,jj)                                                                       &  
     284                       &     + (   ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &  
     285                       &         - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &  
     286                       &     * 0.5 * r1_e12f(ji,jj)  
     287               END DO  
    267288            END DO 
    268289         END SELECT 
    269290 
    270          IF( ln_sco ) THEN 
    271             zwz(:,:) = zwz(:,:) / fse3f(:,:,jk) 
    272             zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 
    273             zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 
    274          ELSE 
    275             zwx(:,:) = e2u(:,:) * un(:,:,jk) 
    276             zwy(:,:) = e1v(:,:) * vn(:,:,jk) 
    277          ENDIF 
     291         zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk)  
     292         zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk)  
     293         zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk)  
    278294 
    279295         ! Compute and add the vorticity term trend 
    … …  
    418434 
    419435 
    420    SUBROUTINE vor_ens( kt, kvor, pua, pva ) 
     436   SUBROUTINE vor_ens( kt, kvor, pun, pvn, pua, pva ) 
    421437      !!---------------------------------------------------------------------- 
    422438      !!                ***  ROUTINE vor_ens  *** 
    … …  
    448464      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend 
    449465      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend 
     466      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities 
    450467      ! 
    451468      INTEGER  ::   ji, jj, jk           ! dummy loop indices 
    … …  
    475492         SELECT CASE( kvor )      ! vorticity considered 
    476493         CASE ( 1 )   ;   zwz(:,:) =                  ff(:,:)      ! planetary vorticity (Coriolis) 
    477          CASE ( 2 )   ;   zwz(:,:) =   rotn(:,:,jk)                ! relative  vorticity 
     494         CASE ( 2 )                                                ! relative  vorticity 
     495            DO jj = 1, jpjm1  
     496               DO ji = 1, fs_jpim1   ! vector opt.  
     497                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &  
     498                     &          - e1u(ji  ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e12f(ji,jj)  
     499               END DO  
     500            END DO 
    478501         CASE ( 3 )                                                ! metric term 
    479             DO jj = 1, jpjm1 
    480                DO ji = 1, fs_jpim1   ! vector opt. 
    481                   zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    482                        &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    483                        &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) 
    484                END DO 
    485             END DO 
    486          CASE ( 4 )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) )    ! total (relative + planetary vorticity) 
     502            DO jj = 1, jpjm1  
     503               DO ji = 1, fs_jpim1   ! vector opt.  
     504                  zwz(ji,jj) = (   ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &  
     505                       &         - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &  
     506                       &     * 0.5 * r1_e12f(ji,jj)  
     507               END DO  
     508            END DO  
     509         CASE ( 4 )                                                ! total (relative + planetary vorticity) 
     510            DO jj = 1, jpjm1  
     511               DO ji = 1, fs_jpim1   ! vector opt.  
     512                  zwz(ji,jj) = ff(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &  
     513                     &                      - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) &  
     514                     &                   * r1_e12f(ji,jj)  
     515               END DO  
     516            END DO 
    487517         CASE ( 5 )                                                ! total (coriolis + metric) 
    488             DO jj = 1, jpjm1 
    489                DO ji = 1, fs_jpim1   ! vector opt. 
    490                   zwz(ji,jj) = ( ff (ji,jj)                                                                       & 
    491                        &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    492                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    493                        &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                                & 
    494                        &       ) 
    495                END DO 
     518            DO jj = 1, jpjm1  
     519               DO ji = 1, fs_jpim1   ! vector opt.  
     520                  zwz(ji,jj) = ff(ji,jj)                                                                         &  
     521                       &     + (   ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &  
     522                       &         - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &  
     523                       &     * 0.5 * r1_e12f(ji,jj)  
     524               END DO  
    496525            END DO 
    497526         END SELECT 
    498          ! 
    499          IF( ln_sco ) THEN 
    500             DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop  
    501                DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking 
    502                   zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk) 
    503                   zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk) 
    504                   zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk) 
    505                END DO 
    506             END DO 
    507          ELSE 
    508             DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop  
    509                DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking 
    510                   zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk) 
    511                   zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk) 
    512                END DO 
    513             END DO 
    514          ENDIF 
    515          ! 
     527         !                                 !==  horizontal fluxes  ==!  
     528         zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk)  
     529         zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk)  
     530         zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk)  
     531         !  
    516532         ! Compute and add the vorticity term trend 
    517533         ! ---------------------------------------- 
    … …  
    536552 
    537553 
    538    SUBROUTINE vor_een( kt, kvor, pua, pva ) 
     554   SUBROUTINE vor_een( kt, kvor, pun, pvn, pua, pva ) 
    539555      !!---------------------------------------------------------------------- 
    540556      !!                ***  ROUTINE vor_een  *** 
    … …  
    559575      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend 
    560576      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend 
     577      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! now velocities 
    561578      !! 
    562579      INTEGER  ::   ji, jj, jk                                    ! dummy loop indices 
    … …  
    604621                     ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   & 
    605622                        &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) ) 
    606                      IF( ze3 /= 0._wp )   ze3f(ji,jj,jk) = 4.0_wp / ze3 
     623                     IF( ze3 /= 0._wp ) THEN  ;   ze3f(ji,jj,jk) = 4.0_wp / ze3 
     624                     ELSE                     ;   ze3f(ji,jj) = 0._wp 
     625                     ENDIF 
    607626                  END DO 
    608627               END DO 
    … …  
    616635                     zmsk = (                   tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   & 
    617636                        &                     + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk) ) 
    618                      IF( ze3 /= 0._wp )   ze3f(ji,jj,jk) = zmsk / ze3 
     637                     IF( ze3 /= 0._wp ) THEN  ;  ze3f(ji,jj,jk) = zmsk / ze3 
     638                     ELSE                     ;  ze3f(ji,jj) = 0._wp 
     639                     ENDIF 
    619640                  END DO 
    620641               END DO 
    … …  
    639660            zwz(:,:) = ff(:,:)      * ze3f(:,:,jk) 
    640661         CASE ( 2 )                                                ! relative  vorticity 
    641             zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk) 
    642          CASE ( 3 )                                                ! metric term 
    643             DO jj = 1, jpjm1 
    644                DO ji = 1, fs_jpim1   ! vector opt. 
    645                   zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    646                        &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    647                        &     * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f(ji,jj,jk) 
    648                END DO 
    649             END DO 
    650             CALL lbc_lnk( zwz, 'F', 1. ) 
     662            DO jj = 1, jpjm1  
     663               DO ji = 1, fs_jpim1   ! vector opt.  
     664                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &  
     665                     &          - e1u(ji  ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) &  
     666                     &       * r1_e12f(ji,jj) * ze3f(ji,jj)  
     667               END DO  
     668            END DO 
     669        CASE ( 3 )                                                ! metric term 
     670            DO jj = 1, jpjm1  
     671               DO ji = 1, fs_jpim1   ! vector opt.  
     672                  zwz(ji,jj) = (   ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &  
     673                       &         - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &  
     674                       &     * 0.5 * r1_e12f(ji,jj) * ze3f(ji,jj)  
     675               END DO  
     676            END DO 
    651677        CASE ( 4 )                                                ! total (relative + planetary vorticity) 
    652             zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk) 
     678            DO jj = 1, jpjm1  
     679               DO ji = 1, fs_jpim1   ! vector opt.  
     680                  zwz(ji,jj) = (  ff(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &  
     681                     &                         - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) &  
     682                     &                      * r1_e12f(ji,jj) ) * ze3f(ji,jj)  
     683               END DO  
     684            END DO 
    653685         CASE ( 5 )                                                ! total (coriolis + metric) 
    654             DO jj = 1, jpjm1 
    655                DO ji = 1, fs_jpim1   ! vector opt. 
    656                   zwz(ji,jj) = ( ff (ji,jj)                                                                       & 
    657                        &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    658                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    659                        &       * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )                                                & 
    660                        &       ) * ze3f(ji,jj,jk) 
    661                END DO 
    662             END DO 
    663             CALL lbc_lnk( zwz, 'F', 1. ) 
     686            DO jj = 1, jpjm1  
     687               DO ji = 1, fs_jpim1   ! vector opt.  
     688                  zwz(ji,jj) = (  ff(ji,jj)                                                                         &  
     689                       &        + (   ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       &  
     690                       &            - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   &  
     691                       &        * 0.5 * r1_e12f(ji,jj)   ) * ze3f(ji,jj)  
     692               END DO  
     693            END DO  
    664694         END SELECT 
    665  
    666          zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 
    667          zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 
     695         ! 
     696         CALL lbc_lnk( zwz, 'F', 1. )  
     697         ! 
     698         !                                   !==  horizontal fluxes  ==! 
     699         zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * pun(:,:,jk) 
     700         zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * pvn(:,:,jk) 
    668701 
    669702         ! Compute and add the vorticity term trend 

  • branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90

    r7470 r7481  
    2424   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2525   USE prtctl          ! Print control 
    26    USE sbcwave,ONLY : cdn_wave !wave module 
    2726 
    2827   IMPLICIT NONE 

  • branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

    r7471 r7481  
    11351135      !                                                      !       Stokes drift u      !  
    11361136      !                                                      ! ========================= !   
    1137          IF( srcv(jpr_sdrftx)%laction ) zusd2dt(:,:) = frcv(jpr_sdrftx)%z3(:,:,1)  
     1137         IF( srcv(jpr_sdrftx)%laction ) ut0sd(:,:) = frcv(jpr_sdrftx)%z3(:,:,1)  
    11381138      !  
    11391139      !                                                      ! ========================= !   
    11401140      !                                                      !       Stokes drift v      !  
    11411141      !                                                      ! ========================= !   
    1142          IF( srcv(jpr_sdrfty)%laction ) zvsd2dt(:,:) = frcv(jpr_sdrfty)%z3(:,:,1)  
     1142         IF( srcv(jpr_sdrfty)%laction ) vt0sd(:,:) = frcv(jpr_sdrfty)%z3(:,:,1)  
    11431143      !  
    11441144      !                                                      ! ========================= !   
    … …  
    11501150      !                                                      !  Significant wave height  !  
    11511151      !                                                      ! ========================= !   
    1152          IF( srcv(jpr_hsig)%laction ) swh(:,:) = frcv(jpr_hsig)%z3(:,:,1)  
     1152         IF( srcv(jpr_hsig)%laction ) hsw(:,:) = frcv(jpr_hsig)%z3(:,:,1)  
    11531153      !  
    11541154      !                                                      ! ========================= !   
    … …  
    11591159         ! Calculate the 3D Stokes drift both in coupled and not fully uncoupled mode  
    11601160         IF( srcv(jpr_sdrftx)%laction .OR. srcv(jpr_sdrfty)%laction .OR. srcv(jpr_wper)%laction &  
    1161                                                                     .OR. srcv(jpr_hsig)%laction ) THEN  
     1161                                                                    .OR. srcv(jpr_hsig)%laction ) &  
    11621162            CALL sbc_stokes()  
    1163             IF( ln_zdfqiao .AND. .NOT. srcv(jpr_wnum)%laction ) CALL sbc_qiao()  
    1164          ENDIF  
    1165          IF( ln_zdfqiao .AND. srcv(jpr_wnum)%laction ) CALL sbc_qiao()  
    11661163      ENDIF  
    11671164      !                                                      ! ========================= !   

  • branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r7471 r7481  
    315315 
    316316      IF( nn_ice == 4      )   CALL cice_sbc_init( nsbc )      ! CICE initialisation 
     317      !  
     318      IF( ln_wave          )   CALL sbc_wave_init              ! surface wave initialisation  
     319      !  
    317320       
    318321   END SUBROUTINE sbc_init 

  • branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90

    r7471 r7481  
    77   !!         :  3.4  !   2012-10  (Adani M)                 Stokes Drift  
    88   !!            3.6  !   2014-09  (Clementi E, Oddo P)New Stokes Drift Computation 
    9    !!---------------------------------------------------------------------- 
    10  
    11    !!---------------------------------------------------------------------- 
     9   !!             -   !   2016-12  (G. Madec, E. Clementi) update Stoke drift computation  
     10   !!                                                     + add sbc_wave_ini routine 
     11   !!---------------------------------------------------------------------- 
     12 
     13   !!---------------------------------------------------------------------- 
     14   !!   sbc_stokes    : calculate 3D Stokes-drift velocities 
    1215   !!   sbc_wave      : wave data from wave model in netcdf files  
    13    !!---------------------------------------------------------------------- 
    14    USE oce            !  
    15    USE sbc_oce       ! Surface boundary condition: ocean fields 
    16    USE bdy_oce        ! 
    17    USE domvvl         ! 
     16   !!   sbc_wave_init : initialisation fo surface waves 
     17   !!---------------------------------------------------------------------- 
     18   USE oce            ! ocean variables 
     19   USE sbc_oce        ! Surface boundary condition: ocean fields 
     20   USE zdf_oce,  ONLY : ln_zdfqiao 
     21   USE bdy_oce        ! open boundary condition variables 
     22   USE domvvl         ! domain: variable volume layers 
    1823   ! 
    1924   USE iom            ! I/O manager library 
    2025   USE in_out_manager ! I/O manager 
    2126   USE lib_mpp        ! distribued memory computing library 
    22    USE fldread       ! read input fields 
     27   USE fldread        ! read input fields 
    2328   USE wrk_nemo       ! 
    2429   USE phycst         ! physical constants  
    … …  
    2732   PRIVATE 
    2833 
    29    PUBLIC   sbc_stokes, sbc_qiao  ! routines called in sbccpl 
    30    PUBLIC   sbc_wave    ! routine called in sbcmod 
     34   PUBLIC   sbc_stokes      ! routine called in sbccpl  
     35   PUBLIC   sbc_wave        ! routine called in sbcmod  
     36   PUBLIC   sbc_wave_init   ! routine called in sbcmod 
    3137    
    3238   ! Variables checking if the wave parameters are coupled (if not, they are read from file) 
    33    LOGICAL, PUBLIC     ::   cpl_hsig=.FALSE. 
    34    LOGICAL, PUBLIC     ::   cpl_phioc=.FALSE. 
    35    LOGICAL, PUBLIC     ::   cpl_sdrftx=.FALSE. 
    36    LOGICAL, PUBLIC     ::   cpl_sdrfty=.FALSE. 
    37    LOGICAL, PUBLIC     ::   cpl_wper=.FALSE. 
    38    LOGICAL, PUBLIC     ::   cpl_wnum=.FALSE. 
    39    LOGICAL, PUBLIC     ::   cpl_wstrf=.FALSE. 
    40    LOGICAL, PUBLIC     ::   cpl_wdrag=.FALSE. 
    41  
    42    INTEGER ::   jpfld                ! number of files to read for stokes drift 
    43    INTEGER ::   jp_usd               ! index of stokes drift  (i-component) (m/s)    at T-point 
    44    INTEGER ::   jp_vsd               ! index of stokes drift  (j-component) (m/s)    at T-point 
    45    INTEGER ::   jp_swh               ! index of significant wave hight      (m)      at T-point 
    46    INTEGER ::   jp_wmp               ! index of mean wave period            (s)      at T-point 
    47  
    48    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_cd    ! structure of input fields (file informations, fields read) Drag Coefficient 
    49    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_sd    ! structure of input fields (file informations, fields read) Stokes Drift 
    50    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_wn    ! structure of input fields (file informations, fields read) wave number for Qiao 
    51    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_tauoc ! structure of input fields (file informations, fields read) normalized wave stress into the ocean 
    52    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: cdn_wave  
    53    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: swh,wmp, wnum 
    54    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: tauoc_wave 
    55    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: tsd2d 
    56    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: zusd2dt, zvsd2dt 
    57    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     :: usd3d, vsd3d, wsd3d  
    58    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     :: usd3dt, vsd3dt 
     39   LOGICAL, PUBLIC ::   cpl_hsig   = .FALSE. 
     40   LOGICAL, PUBLIC ::   cpl_phioc  = .FALSE. 
     41   LOGICAL, PUBLIC ::   cpl_sdrftx = .FALSE. 
     42   LOGICAL, PUBLIC ::   cpl_sdrfty = .FALSE. 
     43   LOGICAL, PUBLIC ::   cpl_wper   = .FALSE. 
     44   LOGICAL, PUBLIC ::   cpl_wnum   = .FALSE. 
     45   LOGICAL, PUBLIC ::   cpl_wstrf  = .FALSE. 
     46   LOGICAL, PUBLIC ::   cpl_wdrag  = .FALSE. 
     47 
     48   INTEGER ::   jpfld    ! number of files to read for stokes drift 
     49   INTEGER ::   jp_usd   ! index of stokes drift  (i-component) (m/s)    at T-point 
     50   INTEGER ::   jp_vsd   ! index of stokes drift  (j-component) (m/s)    at T-point 
     51   INTEGER ::   jp_hsw   ! index of significant wave hight      (m)      at T-point 
     52   INTEGER ::   jp_wmp   ! index of mean wave period            (s)      at T-point 
     53 
     54   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_cd    ! structure of input fields (file informations, fields read) Drag Coefficient 
     55   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_sd    ! structure of input fields (file informations, fields read) Stokes Drift 
     56   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_wn    ! structure of input fields (file informations, fields read) wave number for Qiao 
     57   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_tauoc ! structure of input fields (file informations, fields read) normalized wave stress into the ocean 
     58   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   cdn_wave            !: 
     59   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   hsw, wmp, wnum      !: 
     60   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauoc_wave          !: 
     61   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tsd2d               !: 
     62   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   div_sd              !: barotropic stokes drift divergence 
     63   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   ut0sd, vt0sd        !: surface Stokes drift velocities at t-point 
     64   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) ::   usd  , vsd  , wsd   !: Stokes drift velocities at u-, v- & w-points, resp. 
    5965 
    6066   !! * Substitutions 
    61 #  include "domzgr_substitute.h90" 
    6267#  include "vectopt_loop_substitute.h90" 
    6368   !!---------------------------------------------------------------------- 
    … …  
    8085      !! ** action   
    8186      !!--------------------------------------------------------------------- 
    82       INTEGER                ::   jj,ji,jk  
    83       REAL(wp)                       ::  ztransp, zfac, zsp0, zk, zus, zvs 
    84       REAL(wp), DIMENSION(:,:,:), POINTER :: ze3hdiv   ! 3D workspace 
    85       !!--------------------------------------------------------------------- 
    86       ! 
    87  
    88       CALL wrk_alloc( jpi,jpj,jpk, ze3hdiv ) 
    89       DO jk = 1, jpk 
    90          DO jj = 1, jpj 
    91             DO ji = 1, jpi 
    92                ! On T grid 
    93                ! Stokes transport speed estimated from Hs and Tmean 
    94                ztransp = 2.0_wp*rpi*swh(ji,jj)**2.0_wp/(16.0_wp*MAX(wmp(ji,jj),0.0000001_wp)) 
    95                ! Stokes surface speed 
    96                zsp0 = SQRT( zusd2dt(ji,jj)**2 + zvsd2dt(ji,jj)**2) 
    97                ! Wavenumber scale 
    98                zk = ABS(zsp0)/MAX(ABS(5.97_wp*ztransp),0.0000001_wp) 
    99                ! Depth attenuation 
    100                zfac = EXP(-2.0_wp*zk*fsdept(ji,jj,jk))/(1.0_wp+8.0_wp*zk*fsdept(ji,jj,jk)) 
     87      INTEGER  ::   jj, ji, jk   ! dummy loop argument 
     88      INTEGER  ::   ik           ! local integer  
     89      REAL(wp) ::  ztransp, zfac, ztemp 
     90      REAL(wp) ::  zdep_u, zdep_v, zkh_u, zkh_v, zda_u, zda_v 
     91      REAL(wp), DIMENSION(:,:)  , POINTER ::   zk_t, zk_u, zk_v, zu0_sd, zv0_sd   ! 2D workspace 
     92      REAL(wp), DIMENSION(:,:,:), POINTER ::   ze3divh                            ! 3D workspace 
     93 
     94      !!--------------------------------------------------------------------- 
     95      ! 
     96 
     97      CALL wrk_alloc( jpi,jpj,jpk,   ze3divh ) 
     98      CALL wrk_alloc( jpi,jpj,       zk_t, zk_u, zk_v, zu0_sd, zv0_sd )  
     99      ! 
     100      zfac = 2.0_wp * rpi / 16.0_wp 
     101      DO jj = 1, jpj               ! exp. wave number at t-point    (Eq. (19) in Breivick et al. (2014) ) 
     102         DO ji = 1, jpi 
     103            ! Stokes drift velocity estimated from Hs and Tmean 
     104            ztransp = zfac * hsw(ji,jj)*hsw(ji,jj) / MAX( wmp(ji,jj), 0.0000001_wp ) 
     105            ! Stokes surface speed 
     106            tsd2d(ji,jj) = SQRT( ut0sd(ji,jj)*ut0sd(ji,jj) + vt0sd(ji,jj)*vt0sd(ji,jj)) 
     107            ! Wavenumber scale 
     108            zk_t(ji,jj) = ABS( tsd2d(ji,jj) ) / MAX( ABS( 5.97_wp*ztransp ), 0.0000001_wp ) 
     109         END DO 
     110      END DO 
     111      DO jj = 1, jpjm1              ! exp. wave number & Stokes drift velocity at u- & v-points 
     112         DO ji = 1, jpim1 
     113            zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) ) 
     114            zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) ) 
     115            ! 
     116            zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) ) 
     117            zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) ) 
     118         END DO 
     119      END DO 
     120      ! 
     121      !                       !==  horizontal Stokes Drift 3D velocity  ==! 
     122      DO jk = 1, jpkm1 
     123         DO jj = 2, jpjm1 
     124            DO ji = 2, jpim1 
     125               zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) ) 
     126               zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) ) 
     127               !                           
     128               zkh_u = zk_u(ji,jj) * zdep_u     ! k * depth 
     129               zkh_v = zk_v(ji,jj) * zdep_v 
     130               !                                ! Depth attenuation 
     131               zda_u = EXP( -2.0_wp*zkh_u ) / ( 1.0_wp + 8.0_wp*zkh_u ) 
     132               zda_v = EXP( -2.0_wp*zkh_v ) / ( 1.0_wp + 8.0_wp*zkh_v ) 
    101133               ! 
    102                usd3dt(ji,jj,jk) = zfac * zusd2dt(ji,jj) * tmask(ji,jj,jk) 
    103                vsd3dt(ji,jj,jk) = zfac * zvsd2dt(ji,jj) * tmask(ji,jj,jk) 
     134               usd(ji,jj,jk) = zda_u * zk_u(ji,jj) * umask(ji,jj,jk) 
     135               vsd(ji,jj,jk) = zda_v * zk_v(ji,jj) * vmask(ji,jj,jk) 
    104136            END DO 
    105137         END DO 
    106       END DO  
    107       ! Into the U and V Grid 
    108       DO jk = 1, jpkm1 
    109          DO jj = 1, jpjm1 
    110             DO ji = 1, fs_jpim1 
    111                usd3d(ji,jj,jk) = 0.5 *  umask(ji,jj,jk) *   & 
    112                                &  ( usd3dt(ji,jj,jk) + usd3dt(ji+1,jj,jk) ) 
    113                vsd3d(ji,jj,jk) = 0.5 *  vmask(ji,jj,jk) *   & 
    114                                &  ( vsd3dt(ji,jj,jk) + vsd3dt(ji,jj+1,jk) ) 
    115             END DO 
    116          END DO 
    117       END DO 
    118       ! 
    119       CALL lbc_lnk( usd3d(:,:,:), 'U', -1. ) 
    120       CALL lbc_lnk( vsd3d(:,:,:), 'V', -1. ) 
    121       ! 
    122       DO jk = 1, jpkm1               ! Horizontal divergence 
     138      END DO 
     139      CALL lbc_lnk( usd(:,:,:), 'U', vsd(:,:,:), 'V', -1. ) 
     140      ! 
     141      !                       !==  vertical Stokes Drift 3D velocity  ==! 
     142      ! 
     143      DO jk = 1, jpkm1               ! Horizontal e3*divergence 
    123144         DO jj = 2, jpj 
    124145            DO ji = fs_2, jpi 
    125                ze3hdiv(ji,jj,jk) = (  e2u(ji  ,jj) * usd3d(ji  ,jj,jk)     & 
    126                   &                 - e2u(ji-1,jj) * usd3d(ji-1,jj,jk)     & 
    127                   &                 + e1v(ji,jj  ) * vsd3d(ji,jj  ,jk)     & 
    128                   &                 - e1v(ji,jj-1) * vsd3d(ji,jj-1,jk)   ) * r1_e12t(ji,jj) 
     146               ze3divh(ji,jj,jk) = (  e2u(ji  ,jj) * e3u_n(ji  ,jj,jk) * usd(ji,  jj,jk)    & 
     147                  &                 - e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) * usd(ji-1,jj,jk)    & 
     148                  &                 + e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) * vsd(ji,jj  ,jk)    & 
     149                  &                 - e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vsd(ji,jj-1,jk)  ) * r1_e12t(ji,jj) 
    129150            END DO 
    130151         END DO 
    … …  
    132153      ! 
    133154      IF( .NOT. AGRIF_Root() ) THEN 
    134          IF( nbondi ==  1 .OR. nbondi == 2 )   ze3hdiv(nlci-1,   :  ,:) = 0._wp      ! east 
    135          IF( nbondi == -1 .OR. nbondi == 2 )   ze3hdiv(  2   ,   :  ,:) = 0._wp      ! west 
    136          IF( nbondj ==  1 .OR. nbondj == 2 )   ze3hdiv(  :   ,nlcj-1,:) = 0._wp      ! north 
    137          IF( nbondj == -1 .OR. nbondj == 2 )   ze3hdiv(  :   ,  2   ,:) = 0._wp      ! south 
    138       ENDIF 
    139       ! 
    140       CALL lbc_lnk( ze3hdiv, 'T', 1. ) 
    141       ! 
    142       DO jk = jpkm1, 1, -1                   ! integrate from the bottom the e3t * hor. divergence 
    143          wsd3d(:,:,jk) = wsd3d(:,:,jk+1) - fse3t_n(:,:,jk) * ze3hdiv(:,:,jk) 
     155         IF( nbondi ==  1 .OR. nbondi == 2 )   ze3divh(nlci-1,   :  ,:) = 0._wp      ! east 
     156         IF( nbondi == -1 .OR. nbondi == 2 )   ze3divh(  2   ,   :  ,:) = 0._wp      ! west 
     157         IF( nbondj ==  1 .OR. nbondj == 2 )   ze3divh(  :   ,nlcj-1,:) = 0._wp      ! north 
     158         IF( nbondj == -1 .OR. nbondj == 2 )   ze3divh(  :   ,  2   ,:) = 0._wp      ! south 
     159      ENDIF 
     160      ! 
     161      CALL lbc_lnk( ze3divh, 'T', 1. ) 
     162      ! 
     163      IF( .NOT. lk_vvl ) THEN   ;   ik = 1   ! none zero velocity through the sea surface 
     164      ELSE                      ;   ik = 2   ! w=0 at the surface (set one for all in sbc_wave_init) 
     165      ENDIF 
     166      DO jk = jpkm1, ik, -1          ! integrate from the bottom the hor. divergence (NB: at k=jpk w is always zero) 
     167         wsd(:,:,jk) = wsd(:,:,jk+1) - ze3divh(:,:,jk) 
    144168      END DO 
    145169#if defined key_bdy 
    146170      IF( lk_bdy ) THEN 
    147171         DO jk = 1, jpkm1 
    148             wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:) 
     172            wsd(:,:,jk) = wsd(:,:,jk) * bdytmask(:,:) 
    149173         END DO 
    150174      ENDIF 
    151175#endif 
    152       CALL wrk_dealloc( jpi,jpj,jpk, ze3hdiv ) 
     176      !                       !==  Horizontal divergence of barotropic Stokes transport  ==! 
     177      div_sd(:,:) = 0._wp 
     178      DO jk = 1, jpkm1                                 !  
     179        div_sd(:,:) = div_sd(:,:) + ze3divh(:,:,jk) 
     180      END DO 
     181      ! 
     182      CALL wrk_dealloc( jpi,jpj,jpk,   ze3divh ) 
     183      CALL wrk_dealloc( jpi,jpj,       zk_t, zk_u, zk_v, zu0_sd, zv0_sd ) 
    153184      ! 
    154185   END SUBROUTINE sbc_stokes 
    155186 
    156    SUBROUTINE sbc_qiao 
    157       !!--------------------------------------------------------------------- 
    158       !!                     ***  ROUTINE sbc_qiao  *** 
    159       !! 
    160       !! ** Purpose :   Qiao formulation for wave enhanced turbulence 
    161       !!                2010 (DOI: 10.1007/s10236-010-0326)  
    162       !! 
    163       !! ** Method  : -  
    164       !! ** action   
    165       !!--------------------------------------------------------------------- 
    166       INTEGER :: jj, ji 
    167  
    168       ! Calculate the module of the stokes drift on T grid 
    169       !------------------------------------------------- 
    170       DO jj = 1, jpj 
    171          DO ji = 1, jpi 
    172             tsd2d(ji,jj) = SQRT( zusd2dt(ji,jj) * zusd2dt(ji,jj) + zvsd2dt(ji,jj) * zvsd2dt(ji,jj) ) 
    173          END DO 
    174       END DO 
    175       ! 
    176    END SUBROUTINE sbc_qiao 
    177187 
    178188   SUBROUTINE sbc_wave( kt ) 
    … …  
    190200      !! ** action   
    191201      !!--------------------------------------------------------------------- 
    192       USE zdf_oce,  ONLY : ln_zdfqiao 
    193  
    194       IMPLICIT NONE 
    195  
    196       INTEGER, INTENT( in  ) ::   kt       ! ocean time step 
    197       ! 
    198       INTEGER                ::   ierror   ! return error code 
    199       INTEGER                ::   ifpr 
    200       INTEGER                ::   ios      ! Local integer output status for namelist read 
    201       ! 
     202      INTEGER, INTENT(in   ) ::   kt   ! ocean time step 
     203      !!--------------------------------------------------------------------- 
     204      ! 
     205      IF( ln_cdgw .AND. .NOT. cpl_wdrag ) THEN     !==  Neutral drag coefficient  ==! 
     206         CALL fld_read( kt, nn_fsbc, sf_cd )             ! read from external forcing 
     207         cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1) 
     208      ENDIF 
     209 
     210      IF( ln_tauoc .AND. .NOT. cpl_wstrf ) THEN    !==  Wave induced stress  ==! 
     211         CALL fld_read( kt, nn_fsbc, sf_tauoc )          ! read wave norm stress from external forcing 
     212         tauoc_wave(:,:) = sf_tauoc(1)%fnow(:,:,1) 
     213      ENDIF 
     214 
     215      IF( ln_sdw )  THEN                           !==  Computation of the 3d Stokes Drift  ==!  
     216         ! 
     217         IF( jpfld > 0 ) THEN                            ! Read from file only if the field is not coupled 
     218            CALL fld_read( kt, nn_fsbc, sf_sd )          ! read wave parameters from external forcing 
     219            IF( jp_hsw > 0 )   hsw  (:,:) = sf_sd(jp_hsw)%fnow(:,:,1)   ! significant wave height 
     220            IF( jp_wmp > 0 )   wmp  (:,:) = sf_sd(jp_wmp)%fnow(:,:,1)   ! wave mean period 
     221            IF( jp_usd > 0 )   ut0sd(:,:) = sf_sd(jp_usd)%fnow(:,:,1)   ! 2D zonal Stokes Drift at T point 
     222            IF( jp_vsd > 0 )   vt0sd(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)   ! 2D meridional Stokes Drift at T point 
     223         ENDIF 
     224         ! 
     225         ! Read also wave number if needed, so that it is available in coupling routines 
     226         IF( ln_zdfqiao .AND. .NOT.cpl_wnum ) THEN 
     227            CALL fld_read( kt, nn_fsbc, sf_wn )          ! read wave parameters from external forcing 
     228            wnum(:,:) = sf_wn(1)%fnow(:,:,1) 
     229         ENDIF 
     230            
     231         !                                         !==  Computation of the 3d Stokes Drift  ==!  
     232         ! 
     233         IF( jpfld == 4 )   CALL sbc_stokes()            ! Calculate only if required fields are read 
     234         !                                               ! In coupled wave model-NEMO case the call is done after coupling 
     235         ! 
     236      ENDIF 
     237      ! 
     238   END SUBROUTINE sbc_wave 
     239 
     240 
     241   SUBROUTINE sbc_wave_init 
     242      !!--------------------------------------------------------------------- 
     243      !!                     ***  ROUTINE sbc_wave_init  *** 
     244      !! 
     245      !! ** Purpose :   read wave parameters from wave model  in netcdf files. 
     246      !! 
     247      !! ** Method  : - Read namelist namsbc_wave 
     248      !!              - Read Cd_n10 fields in netcdf files  
     249      !!              - Read stokes drift 2d in netcdf files  
     250      !!              - Read wave number in netcdf files  
     251      !!              - Compute 3d stokes drift using Breivik et al.,2014 
     252      !!                formulation 
     253      !! ** action   
     254      !!--------------------------------------------------------------------- 
     255      INTEGER ::   ierror, ios   ! local integer 
     256      INTEGER ::   ifpr 
     257      !! 
    202258      CHARACTER(len=100)     ::  cn_dir                          ! Root directory for location of drag coefficient files 
    203259      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i     ! array of namelist informations on the fields to read 
    204260      TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd,  & 
    205                              &   sn_swh, sn_wmp, sn_wnum, sn_tauoc      ! informations about the fields to be read 
    206       !! 
    207       NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_swh, sn_wmp, sn_wnum, sn_tauoc 
    208       !!--------------------------------------------------------------------- 
    209       ! 
    210       !                                         ! -------------------- ! 
    211       IF( kt == nit000 ) THEN                   ! First call kt=nit000 ! 
    212          !                                      ! -------------------- ! 
    213          REWIND( numnam_ref )              ! Namelist namsbc_wave in reference namelist : File for drag coeff. from wave model 
    214          READ  ( numnam_ref, namsbc_wave, IOSTAT = ios, ERR = 901) 
    215 901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in reference namelist', lwp ) 
    216  
    217          REWIND( numnam_cfg )              ! Namelist namsbc_wave in configuration namelist : File for drag coeff. from wave model 
    218          READ  ( numnam_cfg, namsbc_wave, IOSTAT = ios, ERR = 902 ) 
    219 902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in configuration namelist', lwp ) 
    220          IF(lwm) WRITE ( numond, namsbc_wave ) 
    221          ! 
    222          IF( ln_cdgw ) THEN 
    223             IF( .NOT. cpl_wdrag ) THEN 
    224                ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg 
    225                IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
    226                ! 
    227                                       ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1)   ) 
    228                IF( sn_cdg%ln_tint )   ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) ) 
    229                CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 
    230             ENDIF 
    231             ALLOCATE( cdn_wave(jpi,jpj) ) 
    232             cdn_wave(:,:) = 0.0 
    233          ENDIF 
    234  
    235          IF( ln_tauoc ) THEN 
    236             IF( .NOT. cpl_wstrf ) THEN 
    237                ALLOCATE( sf_tauoc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_tauoc 
    238                IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
    239                ! 
    240                                        ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1)   ) 
    241                IF( sn_tauoc%ln_tint )  ALLOCATE( sf_tauoc(1)%fdta(jpi,jpj,1,2) ) 
    242                CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' ) 
    243             ENDIF 
    244             ALLOCATE( tauoc_wave(jpi,jpj) ) 
    245          ENDIF 
    246  
    247          IF( ln_sdw ) THEN 
    248             ! Find out how many fields have to be read from file if not coupled 
    249             jpfld=0 
    250             jp_usd=0; jp_vsd=0; jp_swh=0; jp_wmp=0 
    251             IF( .NOT. cpl_sdrftx ) THEN 
    252                jpfld=jpfld+1 
    253                jp_usd=jpfld 
    254             ENDIF 
    255             IF( .NOT. cpl_sdrfty ) THEN 
    256                jpfld=jpfld+1 
    257                jp_vsd=jpfld 
    258             ENDIF 
    259             IF( .NOT. cpl_hsig ) THEN 
    260                jpfld=jpfld+1 
    261                jp_swh=jpfld 
    262             ENDIF 
    263             IF( .NOT. cpl_wper ) THEN 
    264                jpfld=jpfld+1 
    265                jp_wmp=jpfld 
    266             ENDIF 
    267  
    268             ! Read from file only the non-coupled fields  
    269             IF( jpfld > 0 ) THEN 
    270                ALLOCATE( slf_i(jpfld) ) 
    271                IF( jp_usd > 0 ) slf_i(jp_usd) = sn_usd 
    272                IF( jp_vsd > 0 ) slf_i(jp_vsd) = sn_vsd 
    273                IF( jp_swh > 0 ) slf_i(jp_swh) = sn_swh 
    274                IF( jp_wmp > 0 ) slf_i(jp_wmp) = sn_wmp 
    275                ALLOCATE( sf_sd(jpfld), STAT=ierror )           !* allocate and fill sf_sd with stokes drift 
    276                IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
    277                ! 
    278                DO ifpr= 1, jpfld 
    279                   ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) ) 
    280                   IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) ) 
    281                END DO 
    282  
    283                CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 
    284             ENDIF 
    285             ALLOCATE( usd3d(jpi,jpj,jpk),vsd3d(jpi,jpj,jpk),wsd3d(jpi,jpj,jpk) ) 
    286             ALLOCATE( usd3dt(jpi,jpj,jpk),vsd3dt(jpi,jpj,jpk) ) 
    287             ALLOCATE( swh(jpi,jpj), wmp(jpi,jpj) ) 
    288             ALLOCATE( zusd2dt(jpi,jpj), zvsd2dt(jpi,jpj) ) 
    289             usd3d(:,:,:) = 0._wp 
    290             vsd3d(:,:,:) = 0._wp 
    291             wsd3d(:,:,:) = 0._wp 
    292             IF( ln_zdfqiao ) THEN     !==  Vertical mixing enhancement using Qiao,2010  ==! 
    293                IF( .NOT. cpl_wnum ) THEN 
    294                   ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fill sf_wave with sn_wnum 
    295                   IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable toallocate sf_wave structure' ) 
    296                                          ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1)   ) 
    297                   IF( sn_wnum%ln_tint )  ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) ) 
    298                   CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' ) 
    299                ENDIF 
    300                ALLOCATE( wnum(jpi,jpj),tsd2d(jpi,jpj) ) 
    301             ENDIF 
    302          ENDIF 
    303       ENDIF 
    304       ! 
    305       IF( ln_cdgw .AND. .NOT. cpl_wdrag ) THEN              !==  Neutral drag coefficient  ==! 
    306          CALL fld_read( kt, nn_fsbc, sf_cd )      ! read from external forcing 
    307          cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1) 
    308       ENDIF 
    309  
    310       IF( ln_tauoc .AND. .NOT. cpl_wstrf ) THEN             !==  Wave induced stress  ==! 
    311          CALL fld_read( kt, nn_fsbc, sf_tauoc )      !* read wave norm stress from external forcing 
    312          tauoc_wave(:,:) = sf_tauoc(1)%fnow(:,:,1) 
    313       ENDIF 
    314  
    315       IF( ln_sdw )  THEN                         !==  Computation of the 3d Stokes Drift  ==!  
    316          ! 
    317          ! Read from file only if the field is not coupled 
     261                             &   sn_hsw, sn_wmp, sn_wnum, sn_tauoc      ! informations about the fields to be read 
     262      ! 
     263      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wnum, sn_tauoc 
     264      !!--------------------------------------------------------------------- 
     265      ! 
     266      REWIND( numnam_ref )              ! Namelist namsbc_wave in reference namelist : File for drag coeff. from wave model 
     267      READ  ( numnam_ref, namsbc_wave, IOSTAT = ios, ERR = 901) 
     268901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in reference namelist', lwp ) 
     269          
     270      REWIND( numnam_cfg )              ! Namelist namsbc_wave in configuration namelist : File for drag coeff. from wave model 
     271      READ  ( numnam_cfg, namsbc_wave, IOSTAT = ios, ERR = 902 ) 
     272902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in configuration namelist', lwp ) 
     273      IF(lwm) WRITE ( numond, namsbc_wave ) 
     274      ! 
     275      IF( ln_cdgw ) THEN 
     276         IF( .NOT. cpl_wdrag ) THEN 
     277            ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg 
     278            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' ) 
     279            ! 
     280                                   ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1)   ) 
     281            IF( sn_cdg%ln_tint )   ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) ) 
     282            CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave_init', 'Wave module ', 'namsbc_wave' ) 
     283         ENDIF 
     284         ALLOCATE( cdn_wave(jpi,jpj) ) 
     285      ENDIF 
     286 
     287      IF( ln_tauoc ) THEN 
     288         IF( .NOT. cpl_wstrf ) THEN 
     289            ALLOCATE( sf_tauoc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_tauoc 
     290            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' ) 
     291            ! 
     292                                    ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1)   ) 
     293            IF( sn_tauoc%ln_tint )  ALLOCATE( sf_tauoc(1)%fdta(jpi,jpj,1,2) ) 
     294            CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave_init', 'Wave module', 'namsbc_wave' ) 
     295         ENDIF 
     296         ALLOCATE( tauoc_wave(jpi,jpj) ) 
     297      ENDIF 
     298 
     299      IF( ln_sdw ) THEN   ! Find out how many fields have to be read from file if not coupled 
     300         jpfld=0 
     301         jp_usd=0   ;   jp_vsd=0   ;   jp_hsw=0   ;   jp_wmp=0 
     302         IF( .NOT. cpl_sdrftx ) THEN 
     303            jpfld  = jpfld + 1 
     304            jp_usd = jpfld 
     305         ENDIF 
     306         IF( .NOT. cpl_sdrfty ) THEN 
     307            jpfld  = jpfld + 1 
     308            jp_vsd = jpfld 
     309         ENDIF 
     310         IF( .NOT. cpl_hsig ) THEN 
     311            jpfld  = jpfld + 1 
     312            jp_hsw = jpfld 
     313         ENDIF 
     314         IF( .NOT. cpl_wper ) THEN 
     315            jpfld  = jpfld + 1 
     316            jp_wmp = jpfld 
     317         ENDIF 
     318 
     319         ! Read from file only the non-coupled fields  
    318320         IF( jpfld > 0 ) THEN 
    319             CALL fld_read( kt, nn_fsbc, sf_sd )      !* read wave parameters from external forcing 
    320             IF( jp_swh > 0 ) swh(:,:)     = sf_sd(jp_swh)%fnow(:,:,1)   ! significant wave height 
    321             IF( jp_wmp > 0 ) wmp(:,:)     = sf_sd(jp_wmp)%fnow(:,:,1)   ! wave mean period 
    322             IF( jp_usd > 0 ) zusd2dt(:,:) = sf_sd(jp_usd)%fnow(:,:,1)   ! 2D zonal Stokes Drift at T point 
    323             IF( jp_vsd > 0 ) zvsd2dt(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)   ! 2D meridional Stokes Drift at T point 
    324          ENDIF 
    325          ! 
    326          ! Read also wave number if needed, so that it is available in coupling routines 
    327          IF( ln_zdfqiao .AND. .NOT. cpl_wnum ) THEN 
    328             CALL fld_read( kt, nn_fsbc, sf_wn )      !* read wave parameters from external forcing 
    329             wnum(:,:) = sf_wn(1)%fnow(:,:,1) 
    330          ENDIF 
    331             
    332          !==  Computation of the 3d Stokes Drift according to Breivik et al.,2014 
    333          !(DOI: 10.1175/JPO-D-14-0020.1)==!  
    334          ! 
    335          ! Calculate only if no necessary fields are coupled, if not calculate later after coupling 
    336          IF( jpfld == 4 ) THEN 
    337             CALL sbc_stokes() 
    338             IF( ln_zdfqiao .AND. .NOT. cpl_wnum ) THEN 
    339                CALL sbc_qiao() 
    340             ENDIF 
    341          ENDIF 
    342       ENDIF 
    343       ! 
    344    END SUBROUTINE sbc_wave 
    345        
     321            ALLOCATE( slf_i(jpfld) ) 
     322            IF( jp_usd > 0 )   slf_i(jp_usd) = sn_usd 
     323            IF( jp_vsd > 0 )   slf_i(jp_vsd) = sn_vsd 
     324            IF( jp_hsw > 0 )   slf_i(jp_hsw) = sn_hsw 
     325            IF( jp_wmp > 0 )   slf_i(jp_wmp) = sn_wmp 
     326            ALLOCATE( sf_sd(jpfld), STAT=ierror )   !* allocate and fill sf_sd with stokes drift 
     327            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' ) 
     328            ! 
     329            DO ifpr= 1, jpfld 
     330               ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) ) 
     331               IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) ) 
     332            END DO 
     333            ! 
     334            CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave_init', 'Wave module ', 'namsbc_wave' ) 
     335         ENDIF 
     336         ALLOCATE( usd  (jpi,jpj,jpk), vsd  (jpi,jpj,jpk), wsd(jpi,jpj,jpk) ) 
     337         ALLOCATE( hsw  (jpi,jpj)    , wmp  (jpi,jpj)     ) 
     338         ALLOCATE( ut0sd(jpi,jpj)    , vt0sd(jpi,jpj)     ) 
     339         ALLOCATE( div_sd(jpi,jpj) ) 
     340         ALLOCATE( tsd2d (jpi,jpj) ) 
     341         usd(:,:,:) = 0._wp 
     342         vsd(:,:,:) = 0._wp 
     343         wsd(:,:,:) = 0._wp 
     344         ! Wave number needed only if ln_zdfqiao=T 
     345         IF( .NOT. cpl_wnum ) THEN 
     346            ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fill sf_wave with sn_wnum 
     347            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable toallocate sf_wave structure' ) 
     348                                   ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1)   ) 
     349            IF( sn_wnum%ln_tint )  ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) ) 
     350            CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' ) 
     351         ENDIF 
     352         ALLOCATE( wnum(jpi,jpj) ) 
     353      ENDIF 
     354      ! 
     355   END SUBROUTINE sbc_wave_init 
     356 
    346357   !!====================================================================== 
    347358END MODULE sbcwave 

  • branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90

    r7471 r7481  
    101101      !                                               !==  effective transport  ==! 
    102102      IF(ln_wave .AND. ln_sdw) THEN  
    103          DO jk = 1, jpkm1  
    104             zun(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) *      &  
    105                         &  ( un(:,:,jk) + usd3d(:,:,jk) )                     !eulerian transport + Stokes Drift  
    106             zvn(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) *      &  
    107                         &  ( vn(:,:,jk) + vsd3d(:,:,jk) )  
    108             zwn(:,:,jk) = e1e2t(:,:) *                    &  
    109                         &  ( wn(:,:,jk) + wsd3d(:,:,jk) )  
     103         DO jk = 1, jpkm1                                                                      ! eulerian transport + Stokes Drift  
     104            zun(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * ( un(:,:,jk) + usd(:,:,jk) ) 
     105            zvn(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * ( vn(:,:,jk) + vsd(:,:,jk) )  
     106            zwn(:,:,jk) = e1e2t(:,:)               * ( wn(:,:,jk) + wsd(:,:,jk) )  
    110107         END DO  
    111108      ELSE  

  • branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90

    r7470 r7481  
    2424   USE phycst         ! physical constants 
    2525   USE zdfmxl         ! mixed layer 
     26   USE sbcwave ,  ONLY: hsw   ! significant wave height  
     27   !  
    2628   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    2729   USE lib_mpp        ! MPP manager 
    … …  
    197199         zdep(:,:)  = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(:,:),rsmall))))             ! Wave age (eq. 10) 
    198200         zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 
    199       ! 
     201      CASE ( 3 )             ! Roughness given by the wave model (coupled or read in file)  
     202         zhsro(:,:) = hsw(:,:) 
    200203      END SELECT 
    201204 
    … …  
    909912 
    910913      !                                !* Check of some namelist values 
    911       IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' ) 
    912       IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' ) 
    913       IF( nn_z0_met < 0 .OR. nn_z0_met > 2 ) CALL ctl_stop( 'bad flag: nn_z0_met is 0, 1 or 2' ) 
    914       IF( nn_stab_func  < 0 .OR. nn_stab_func  > 3 ) CALL ctl_stop( 'bad flag: nn_stab_func is 0, 1, 2 and 3' ) 
    915       IF( nn_clos       < 0 .OR. nn_clos       > 3 ) CALL ctl_stop( 'bad flag: nn_clos is 0, 1, 2 or 3' ) 
     914      IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' )  
     915      IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' )  
     916      IF( nn_z0_met < 0 .OR. nn_z0_met > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_z0_met is 0, 1, 2 or 3' )  
     917      IF( nn_z0_met == 3 .AND. .NOT.ln_wave ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_wave=T' )  
     918      IF( nn_stab_func  < 0 .OR. nn_stab_func  > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_stab_func is 0, 1, 2 and 3' )  
     919      IF( nn_clos       < 0 .OR. nn_clos       > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' ) 
    916920 
    917921      SELECT CASE ( nn_clos )          !* set the parameters for the chosen closure 

  • branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90

    r7471 r7481  
    7171         DO jj = 1, jpjm1 
    7272            DO ji = 1, fs_jpim1 
    73                qbv(ji,jj,jk) = 1.0 * 0.353553 * swh(ji,jj) * tsd2d(ji,jj) *           & 
     73               qbv(ji,jj,jk) = 1.0 * 0.353553 * hsw(ji,jj) * tsd2d(ji,jj) *           & 
    7474            &                  EXP(3.0 * wnum(ji,jj) *                                &                      
    7575            &                  (-MIN( fsdepw(ji  ,jj  ,jk), fsdepw(ji+1,jj  ,jk),     & 

  • branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/step.F90

    r7471 r7481  
    212212                                  CALL dyn_adv      ( kstp )   ! advection (vector or flux form) 
    213213                                  CALL dyn_vor      ( kstp )   ! vorticity term including Coriolis 
    214           IF( ln_wave .AND. ln_sdw .AND. ln_stcor )          &  
    215                         &         CALL dyn_stcor     ( kstp )  ! Stokes-Coriolis forcing 
    216214                                  CALL dyn_ldf      ( kstp )   ! lateral mixing 
    217215          IF( ln_neptsimp )       CALL dyn_nept_cor ( kstp )   ! add Neptune velocities (simplified) 

  • branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

    r7471 r7481  
    5252   USE dynspg_oce       ! surface pressure gradient        (dyn_spg routine) 
    5353   USE dynspg           ! surface pressure gradient        (dyn_spg routine) 
    54    USE dynstcor         ! simp. form of Stokes-Coriolis 
    5554   USE dynnept          ! simp. form of Neptune effect(dyn_nept_cor routine) 
    5655 
Note: See TracChangeset for help on using the changeset viewer.