Changeset 7422


Ignore:
Timestamp:
2016-12-01T18:17:41+01:00 (4 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
Files:
1 deleted
12 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_INGV_UKMO_2016/NEMOGCM/CONFIG/SHARED/namelist_ref

    r7351 r7422  
    552552/ 
    553553!----------------------------------------------------------------------- 
    554 &namsbc_wave   ! External fields from wave model 
     554&namsbc_wave   ! External fields from wave model                        (ln_wave=T) 
    555555!----------------------------------------------------------------------- 
    556556!              !  file name  ! frequency (hours) ! variable     ! time interp. !  clim   ! 'yearly'/ ! weights  ! rotation ! land/sea mask ! 
     
    10461046   rn_hsro       =  0.02   !  Minimum surface roughness 
    10471047   rn_frac_hs    =   1.3   !  Fraction of wave height as roughness (if nn_z0_met=2) 
    1048    nn_z0_met     =     2   !  Method for surface roughness computation (0/1/2) 
     1048   nn_z0_met     =     2   !  Method for surface roughness computation (0/1/2/3) 
     1049   !                             ! =3 requires ln_wave=T 
    10491050   nn_bc_surf    =     1   !  surface condition (0/1=Dir/Neum) 
    10501051   nn_bc_bot     =     1   !  bottom condition (0/1=Dir/Neum) 
  • branches/2016/dev_INGV_UKMO_2016/NEMOGCM/CONFIG/cfg.txt

    r7351 r7422  
    99ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC 
    1010ORCA2_OFF_PISCES OPA_SRC OFF_SRC TOP_SRC 
     11ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 
    1112GYRE OPA_SRC 
    12 ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 
  • 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  =! 
  • branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

    r7359 r7422  
    11301130      !                                                      !       Stokes drift u      ! 
    11311131      !                                                      ! ========================= !  
    1132          IF( srcv(jpr_sdrftx)%laction ) zusd2dt(:,:) = frcv(jpr_sdrftx)%z3(:,:,1) 
     1132         IF( srcv(jpr_sdrftx)%laction ) ut0sd(:,:) = frcv(jpr_sdrftx)%z3(:,:,1) 
    11331133      ! 
    11341134      !                                                      ! ========================= !  
    11351135      !                                                      !       Stokes drift v      ! 
    11361136      !                                                      ! ========================= !  
    1137          IF( srcv(jpr_sdrfty)%laction ) zvsd2dt(:,:) = frcv(jpr_sdrfty)%z3(:,:,1) 
     1137         IF( srcv(jpr_sdrfty)%laction ) vt0sd(:,:) = frcv(jpr_sdrfty)%z3(:,:,1) 
    11381138      ! 
    11391139      !                                                      ! ========================= !  
     
    11451145      !                                                      !  Significant wave height  ! 
    11461146      !                                                      ! ========================= !  
    1147          IF( srcv(jpr_hsig)%laction ) swh(:,:) = frcv(jpr_hsig)%z3(:,:,1) 
     1147         IF( srcv(jpr_hsig)%laction ) hsw(:,:) = frcv(jpr_hsig)%z3(:,:,1) 
    11481148      ! 
    11491149      !                                                      ! ========================= !  
     
    11561156                                                                    .OR. srcv(jpr_hsig)%laction ) THEN 
    11571157            CALL sbc_stokes() 
    1158             IF( ln_zdfqiao .AND. .NOT. srcv(jpr_wnum)%laction ) CALL sbc_qiao() 
    1159          ENDIF 
    1160          IF( ln_zdfqiao .AND. srcv(jpr_wnum)%laction ) CALL sbc_qiao() 
     1158         ENDIF 
    11611159      ENDIF 
    11621160      !                                                      ! ========================= !  
  • branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r7376 r7422  
    313313      IF( nn_ice == 4 )   CALL cice_sbc_init( nsbc )      ! CICE initialisation 
    314314      ! 
     315      IF( ln_wave     )   CALL sbc_wave_init              ! surface wave initialisation 
     316      ! 
    315317   END SUBROUTINE sbc_init 
    316318 
  • branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90

    r7359 r7422  
    44   !! Wave module  
    55   !!====================================================================== 
    6    !! History :  3.3  !   2011-09  (M. Adani)  Original code: Drag Coefficient  
    7    !!         :  3.4  !   2012-10  (M. Adani)  Stokes Drift  
    8    !!            3.6  !   2014-09  (E. Clementi,P. Oddo) New Stokes Drift Computation 
    9    !!---------------------------------------------------------------------- 
    10  
    11    !!---------------------------------------------------------------------- 
     6   !! History :  3.3  !  2011-09  (M. Adani)  Original code: Drag Coefficient  
     7   !!         :  3.4  !  2012-10  (M. Adani)  Stokes Drift  
     8   !!            3.6  !  2014-09  (E. Clementi,P. Oddo) New Stokes Drift Computation 
     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 phycst         ! physical constants  
     19   USE oce            ! ocean variables 
     20   USE sbc_oce        ! Surface boundary condition: ocean fields 
     21   USE zdf_oce,  ONLY : ln_zdfqiao 
     22   USE bdy_oce        ! open boundary condition variables 
     23   USE domvvl         ! domain: variable volume layers 
     24   ! 
    1825   USE iom            ! I/O manager library 
    1926   USE in_out_manager ! I/O manager 
    2027   USE lib_mpp        ! distribued memory computing library 
    21    USE fldread       ! read input fields 
     28   USE fldread        ! read input fields 
    2229   USE wrk_nemo       ! 
    23    USE phycst         ! physical constants  
    2430 
    2531   IMPLICIT NONE 
    2632   PRIVATE 
    2733 
    28    PUBLIC   sbc_stokes, sbc_qiao  ! routines called in sbccpl 
    29    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 
    3037    
    3138   ! Variables checking if the wave parameters are coupled (if not, they are read from file) 
    32    LOGICAL, PUBLIC     ::   cpl_hsig=.FALSE. 
    33    LOGICAL, PUBLIC     ::   cpl_phioc=.FALSE. 
    34    LOGICAL, PUBLIC     ::   cpl_sdrftx=.FALSE. 
    35    LOGICAL, PUBLIC     ::   cpl_sdrfty=.FALSE. 
    36    LOGICAL, PUBLIC     ::   cpl_wper=.FALSE. 
    37    LOGICAL, PUBLIC     ::   cpl_wnum=.FALSE. 
    38    LOGICAL, PUBLIC     ::   cpl_wstrf=.FALSE. 
    39    LOGICAL, PUBLIC     ::   cpl_wdrag=.FALSE. 
    40  
    41    INTEGER ::   jpfld                ! number of files to read for stokes drift 
    42    INTEGER ::   jp_usd               ! index of stokes drift  (i-component) (m/s)    at T-point 
    43    INTEGER ::   jp_vsd               ! index of stokes drift  (j-component) (m/s)    at T-point 
    44    INTEGER ::   jp_swh               ! index of significant wave hight      (m)      at T-point 
    45    INTEGER ::   jp_wmp               ! index of mean wave period            (s)      at T-point 
    46  
    47    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_cd    ! structure of input fields (file informations, fields read) Drag Coefficient 
    48    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_sd    ! structure of input fields (file informations, fields read) Stokes Drift 
    49    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_wn    ! structure of input fields (file informations, fields read) wave number for Qiao 
    50    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_tauoc ! structure of input fields (file informations, fields read) normalized wave stress into the ocean 
    51    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: cdn_wave  
    52    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: swh,wmp, wnum 
    53    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: tauoc_wave 
    54    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: tsd2d 
    55    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)       :: zusd2dt, zvsd2dt 
    56    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)     :: usd3d, vsd3d, wsd3d  
    57    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. 
    5865 
    5966   !! * Substitutions 
     
    7885      !! ** action   
    7986      !!--------------------------------------------------------------------- 
    80       INTEGER                ::   jj,ji,jk  
    81       REAL(wp)                       ::  ztransp, zfac, zsp0, zk, zus, zvs 
    82       REAL(wp), DIMENSION(:,:,:), POINTER :: ze3hdiv   ! 3D workspace 
    83       !!--------------------------------------------------------------------- 
    84       ! 
    85  
    86       CALL wrk_alloc( jpi,jpj,jpk, ze3hdiv ) 
    87       DO jk = 1, jpk 
    88          DO jj = 1, jpj 
    89             DO ji = 1, jpi 
    90                ! On T grid 
    91                ! Stokes transport speed estimated from Hs and Tmean 
    92                ztransp = 2.0_wp*rpi*swh(ji,jj)**2.0_wp/(16.0_wp*MAX(wmp(ji,jj),0.0000001_wp)) 
     87      INTEGER  ::   jj, ji, jk   ! dummy loop argument 
     88      INTEGER  ::   ik           ! local integer  
     89      REAL(wp) ::  ztransp, zfac, ztemp, zsp0 
     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      CALL wrk_alloc( jpi,jpj,jpk,   ze3divh ) 
     96      CALL wrk_alloc( jpi,jpj,       zk_t, zk_u, zk_v, zu0_sd, zv0_sd ) 
     97      ! 
     98      ! 
     99      zfac =  2.0_wp * rpi / 16.0_wp 
     100      DO jj = 1, jpj                ! exp. wave number at t-point    (Eq. (19) in Breivick et al. (2014) ) 
     101         DO ji = 1, jpi 
     102               ! Stokes drift velocity estimated from Hs and Tmean 
     103               ztransp = zfac * hsw(ji,jj)*hsw(ji,jj) / MAX( wmp(ji,jj) , 0.0000001_wp ) 
    93104               ! Stokes surface speed 
    94                zsp0 = SQRT( zusd2dt(ji,jj)**2 + zvsd2dt(ji,jj)**2) 
     105               zsp0 = SQRT( ut0sd(ji,jj)*ut0sd(ji,jj) + vt0sd(ji,jj)*vt0sd(ji,jj) ) 
     106               tsd2d(ji,jj) = zsp0 
    95107               ! Wavenumber scale 
    96                zk = ABS(zsp0)/MAX(ABS(5.97_wp*ztransp),0.0000001_wp) 
    97                ! Depth attenuation 
    98                zfac = EXP(-2.0_wp*zk*gdept_n(ji,jj,jk))/(1.0_wp+8.0_wp*zk*gdept_n(ji,jj,jk)) 
     108               zk_t(ji,jj) = ABS( zsp0 ) / 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 ) 
    99133               ! 
    100                usd3dt(ji,jj,jk) = zfac * zusd2dt(ji,jj) * tmask(ji,jj,jk) 
    101                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) 
    102136            END DO 
    103137         END DO 
    104       END DO  
    105       ! Into the U and V Grid 
    106       DO jk = 1, jpkm1 
    107          DO jj = 1, jpjm1 
    108             DO ji = 1, fs_jpim1 
    109                usd3d(ji,jj,jk) = 0.5 *  umask(ji,jj,jk) *   & 
    110                                &  ( usd3dt(ji,jj,jk) + usd3dt(ji+1,jj,jk) ) 
    111                vsd3d(ji,jj,jk) = 0.5 *  vmask(ji,jj,jk) *   & 
    112                                &  ( vsd3dt(ji,jj,jk) + vsd3dt(ji,jj+1,jk) ) 
    113             END DO 
    114          END DO 
    115       END DO 
    116       ! 
    117       CALL lbc_lnk( usd3d(:,:,:), 'U', -1. ) 
    118       CALL lbc_lnk( vsd3d(:,:,:), 'V', -1. ) 
    119       ! 
    120       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 
    121144         DO jj = 2, jpj 
    122145            DO ji = fs_2, jpi 
    123                ze3hdiv(ji,jj,jk) = (  e2u(ji  ,jj) * usd3d(ji  ,jj,jk)     & 
    124                   &                 - e2u(ji-1,jj) * usd3d(ji-1,jj,jk)     & 
    125                   &                 + e1v(ji,jj  ) * vsd3d(ji,jj  ,jk)     & 
    126                   &                 - e1v(ji,jj-1) * vsd3d(ji,jj-1,jk)   ) * r1_e1e2t(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_e1e2t(ji,jj) 
    127150            END DO 
    128151         END DO 
     
    130153      ! 
    131154      IF( .NOT. AGRIF_Root() ) THEN 
    132          IF( nbondi ==  1 .OR. nbondi == 2 )   ze3hdiv(nlci-1,   :  ,:) = 0._wp      ! east 
    133          IF( nbondi == -1 .OR. nbondi == 2 )   ze3hdiv(  2   ,   :  ,:) = 0._wp      ! west 
    134          IF( nbondj ==  1 .OR. nbondj == 2 )   ze3hdiv(  :   ,nlcj-1,:) = 0._wp      ! north 
    135          IF( nbondj == -1 .OR. nbondj == 2 )   ze3hdiv(  :   ,  2   ,:) = 0._wp      ! south 
    136       ENDIF 
    137       ! 
    138       CALL lbc_lnk( ze3hdiv, 'T', 1. ) 
    139       ! 
    140       DO jk = jpkm1, 1, -1                   ! integrate from the bottom the e3t * hor. divergence 
    141          wsd3d(:,:,jk) = wsd3d(:,:,jk+1) - e3t_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( ln_linssh ) 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 e3t * hor. divergence (NB: at k=jpk w is always zero) 
     167         wsd(:,:,jk) = wsd(:,:,jk+1) - ze3divh(:,:,jk) 
    142168      END DO 
    143169#if defined key_bdy 
    144170      IF( lk_bdy ) THEN 
    145171         DO jk = 1, jpkm1 
    146             wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:) 
     172            wsd(:,:,jk) = wsd(:,:,jk) * bdytmask(:,:) 
    147173         END DO 
    148174      ENDIF 
    149175#endif 
    150       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 ) 
    151184      ! 
    152185   END SUBROUTINE sbc_stokes 
    153186 
    154    SUBROUTINE sbc_qiao 
    155       !!--------------------------------------------------------------------- 
    156       !!                     ***  ROUTINE sbc_qiao  *** 
    157       !! 
    158       !! ** Purpose :   Qiao formulation for wave enhanced turbulence 
    159       !!                2010 (DOI: 10.1007/s10236-010-0326)  
    160       !! 
    161       !! ** Method  : -  
    162       !! ** action   
    163       !!--------------------------------------------------------------------- 
    164       INTEGER :: jj, ji 
    165  
    166       ! Calculate the module of the stokes drift on T grid 
    167       !------------------------------------------------- 
    168       DO jj = 1, jpj 
    169          DO ji = 1, jpi 
    170             tsd2d(ji,jj) = SQRT( zusd2dt(ji,jj) * zusd2dt(ji,jj) + zvsd2dt(ji,jj) * zvsd2dt(ji,jj) ) 
    171          END DO 
    172       END DO 
    173       ! 
    174    END SUBROUTINE sbc_qiao 
    175187 
    176188   SUBROUTINE sbc_wave( kt ) 
     
    188200      !! ** action   
    189201      !!--------------------------------------------------------------------- 
    190       USE zdf_oce,  ONLY : ln_zdfqiao 
    191  
    192       INTEGER, INTENT( in  ) :: kt       ! ocean time step 
    193       ! 
    194       INTEGER                ::   ierror   ! return error code 
    195       INTEGER                ::   ifpr 
    196       INTEGER                ::   ios      ! Local integer output status for namelist read 
    197       ! 
     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      !! 
    198258      CHARACTER(len=100)     ::  cn_dir                          ! Root directory for location of drag coefficient files 
    199259      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i     ! array of namelist informations on the fields to read 
    200260      TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd,  & 
    201                              &   sn_swh, sn_wmp, sn_wnum, sn_tauoc      ! informations about the fields to be read 
    202       !! 
    203       NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_swh, sn_wmp, sn_wnum, sn_tauoc 
    204       !!--------------------------------------------------------------------- 
    205       ! 
    206       !                                         ! -------------------- ! 
    207       IF( kt == nit000 ) THEN                   ! First call kt=nit000 ! 
    208          !                                      ! -------------------- ! 
    209          REWIND( numnam_ref )              ! Namelist namsbc_wave in reference namelist : File for drag coeff. from wave model 
    210          READ  ( numnam_ref, namsbc_wave, IOSTAT = ios, ERR = 901) 
    211 901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in reference namelist', lwp ) 
     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 ) 
    212269          
    213          REWIND( numnam_cfg )              ! Namelist namsbc_wave in configuration namelist : File for drag coeff. from wave model 
    214          READ  ( numnam_cfg, namsbc_wave, IOSTAT = ios, ERR = 902 ) 
    215 902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in configuration namelist', lwp ) 
    216          IF(lwm) WRITE ( numond, namsbc_wave ) 
    217          ! 
    218          IF( ln_cdgw ) THEN 
    219             IF( .NOT. cpl_wdrag ) THEN 
    220                ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg 
    221                IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
    222                ! 
    223                                       ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1)   ) 
    224                IF( sn_cdg%ln_tint )   ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) ) 
    225                CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 
     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: 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', '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: 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', '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  
     320         IF( jpfld > 0 ) THEN 
     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: 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', 'Wave module ', 'namsbc_wave' ) 
     335         ENDIF 
     336         ALLOCATE( usd  (jpi,jpj,jpk), vsd  (jpi,jpj,jpk), wsd(jpi,jpj,jpk) ) 
     337         ALLOCATE( usd  (jpi,jpj,jpk), vsd  (jpi,jpj,jpk) ) 
     338         ALLOCATE( hsw  (jpi,jpj)    , wmp  (jpi,jpj)     ) 
     339         ALLOCATE( ut0sd(jpi,jpj)    , vt0sd(jpi,jpj)     ) 
     340         ALLOCATE( div_sd(jpi,jpj) ) 
     341         usd(:,:,:) = 0._wp 
     342         vsd(:,:,:) = 0._wp 
     343         wsd(:,:,:) = 0._wp 
     344         IF( ln_zdfqiao ) THEN     !==  Vertical mixing enhancement using Qiao,2010  ==! 
     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: 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' ) 
    226351            ENDIF 
    227             ALLOCATE( cdn_wave(jpi,jpj) ) 
    228          ENDIF 
    229  
    230          IF( ln_tauoc ) THEN 
    231             IF( .NOT. cpl_wstrf ) THEN 
    232                ALLOCATE( sf_tauoc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_tauoc 
    233                IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
    234                ! 
    235                                        ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1)   ) 
    236                IF( sn_tauoc%ln_tint )  ALLOCATE( sf_tauoc(1)%fdta(jpi,jpj,1,2) ) 
    237                CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' ) 
    238             ENDIF 
    239             ALLOCATE( tauoc_wave(jpi,jpj) ) 
    240          ENDIF 
    241  
    242          IF( ln_sdw ) THEN 
    243             ! Find out how many fields have to be read from file if not coupled 
    244             jpfld=0 
    245             jp_usd=0; jp_vsd=0; jp_swh=0; jp_wmp=0 
    246             IF( .NOT. cpl_sdrftx ) THEN 
    247                jpfld=jpfld+1 
    248                jp_usd=jpfld 
    249             ENDIF 
    250             IF( .NOT. cpl_sdrfty ) THEN 
    251                jpfld=jpfld+1 
    252                jp_vsd=jpfld 
    253             ENDIF 
    254             IF( .NOT. cpl_hsig ) THEN 
    255                jpfld=jpfld+1 
    256                jp_swh=jpfld 
    257             ENDIF 
    258             IF( .NOT. cpl_wper ) THEN 
    259                jpfld=jpfld+1 
    260                jp_wmp=jpfld 
    261             ENDIF 
    262  
    263             ! Read from file only the non-coupled fields  
    264             IF( jpfld > 0 ) THEN 
    265                ALLOCATE( slf_i(jpfld) ) 
    266                IF( jp_usd > 0 ) slf_i(jp_usd) = sn_usd 
    267                IF( jp_vsd > 0 ) slf_i(jp_vsd) = sn_vsd 
    268                IF( jp_swh > 0 ) slf_i(jp_swh) = sn_swh 
    269                IF( jp_wmp > 0 ) slf_i(jp_wmp) = sn_wmp 
    270                ALLOCATE( sf_sd(jpfld), STAT=ierror )           !* allocate and fill sf_sd with stokes drift 
    271                IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
    272                ! 
    273                DO ifpr= 1, jpfld 
    274                   ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) ) 
    275                   IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) ) 
    276                END DO 
    277  
    278                CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 
    279             ENDIF 
    280             ALLOCATE( usd3d(jpi,jpj,jpk),vsd3d(jpi,jpj,jpk),wsd3d(jpi,jpj,jpk) ) 
    281             ALLOCATE( usd3dt(jpi,jpj,jpk),vsd3dt(jpi,jpj,jpk) ) 
    282             ALLOCATE( swh(jpi,jpj), wmp(jpi,jpj) ) 
    283             ALLOCATE( zusd2dt(jpi,jpj), zvsd2dt(jpi,jpj) ) 
    284             usd3d(:,:,:) = 0._wp 
    285             vsd3d(:,:,:) = 0._wp 
    286             wsd3d(:,:,:) = 0._wp 
    287             IF( ln_zdfqiao ) THEN     !==  Vertical mixing enhancement using Qiao,2010  ==! 
    288                IF( .NOT. cpl_wnum ) THEN 
    289                   ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fill sf_wave with sn_wnum 
    290                   IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable toallocate sf_wave structure' ) 
    291                                          ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1)   ) 
    292                   IF( sn_wnum%ln_tint )  ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) ) 
    293                   CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' ) 
    294                ENDIF 
    295                ALLOCATE( wnum(jpi,jpj),tsd2d(jpi,jpj) ) 
    296             ENDIF 
    297          ENDIF 
    298       ENDIF 
    299       ! 
    300       IF( ln_cdgw .AND. .NOT. cpl_wdrag ) THEN              !==  Neutral drag coefficient  ==! 
    301          CALL fld_read( kt, nn_fsbc, sf_cd )      ! read from external forcing 
    302          cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1) 
    303       ENDIF 
    304  
    305       IF( ln_tauoc .AND. .NOT. cpl_wstrf ) THEN             !==  Wave induced stress  ==! 
    306          CALL fld_read( kt, nn_fsbc, sf_tauoc )      !* read wave norm stress from external forcing 
    307          tauoc_wave(:,:) = sf_tauoc(1)%fnow(:,:,1) 
    308       ENDIF 
    309  
    310       IF( ln_sdw )  THEN                         !==  Computation of the 3d Stokes Drift  ==!  
    311          ! 
    312          ! Read from file only if the field is not coupled 
    313          IF( jpfld > 0 ) THEN 
    314             CALL fld_read( kt, nn_fsbc, sf_sd )      !* read wave parameters from external forcing 
    315             IF( jp_swh > 0 ) swh(:,:)     = sf_sd(jp_swh)%fnow(:,:,1)   ! significant wave height 
    316             IF( jp_wmp > 0 ) wmp(:,:)     = sf_sd(jp_wmp)%fnow(:,:,1)   ! wave mean period 
    317             IF( jp_usd > 0 ) zusd2dt(:,:) = sf_sd(jp_usd)%fnow(:,:,1)   ! 2D zonal Stokes Drift at T point 
    318             IF( jp_vsd > 0 ) zvsd2dt(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)   ! 2D meridional Stokes Drift at T point 
    319          ENDIF 
    320          ! 
    321          ! Read also wave number if needed, so that it is available in coupling routines 
    322          IF( ln_zdfqiao .AND. .NOT. cpl_wnum ) THEN 
    323             CALL fld_read( kt, nn_fsbc, sf_wn )      !* read wave parameters from external forcing 
    324             wnum(:,:) = sf_wn(1)%fnow(:,:,1) 
    325          ENDIF 
    326             
    327          !==  Computation of the 3d Stokes Drift according to Breivik et al.,2014 
    328          !(DOI: 10.1175/JPO-D-14-0020.1)==!  
    329          ! 
    330          ! Calculate only if no necessary fields are coupled, if not calculate later after coupling 
    331          IF( jpfld == 4 ) THEN 
    332             CALL sbc_stokes() 
    333             IF( ln_zdfqiao .AND. .NOT. cpl_wnum ) THEN 
    334                CALL sbc_qiao() 
    335             ENDIF 
    336          ENDIF 
    337       ENDIF 
    338       ! 
    339    END SUBROUTINE sbc_wave 
    340        
     352            ALLOCATE( wnum(jpi,jpj),tsd2d(jpi,jpj) ) 
     353         ENDIF 
     354      ENDIF 
     355      ! 
     356   END SUBROUTINE sbc_wave_init 
     357 
    341358   !!====================================================================== 
    342359END MODULE sbcwave 
  • branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90

    r7359 r7422  
    107107      !                                         !==  effective transport  ==! 
    108108      IF( ln_wave .AND. ln_sdw )  THEN 
    109          DO jk = 1, jpkm1 
    110             zun(:,:,jk) = e2u(:,:) * e3u_n(:,:,jk) *      & 
    111                         &  ( un(:,:,jk) + usd3d(:,:,jk) )                       ! eulerian transport + Stokes Drift 
    112             zvn(:,:,jk) = e1v(:,:) * e3v_n(:,:,jk) *      & 
    113                         &  ( vn(:,:,jk) + vsd3d(:,:,jk) ) 
    114             zwn(:,:,jk) = e1e2t(:,:) *                    & 
    115                         &  ( wn(:,:,jk) + wsd3d(:,:,jk) ) 
     109         DO jk = 1, jpkm1                                                       ! eulerian transport + Stokes Drift 
     110            zun(:,:,jk) = e2u  (:,:) * e3u_n(:,:,jk) * ( un(:,:,jk) + usd(:,:,jk) ) 
     111            zvn(:,:,jk) = e1v  (:,:) * e3v_n(:,:,jk) * ( vn(:,:,jk) + vsd(:,:,jk) ) 
     112            zwn(:,:,jk) = e1e2t(:,:)                 * ( wn(:,:,jk) + wsd(:,:,jk) ) 
    116113         END DO 
    117114      ELSE 
  • branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90

    r7351 r7422  
    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 
     
    194196         zdep(:,:)  = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(:,:),rsmall))))             ! Wave age (eq. 10) 
    195197         zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 
    196       ! 
     198      CASE ( 3 )             ! Roughness given by the wave model (coupled or read in file) 
     199         zhsro(:,:) = hsw(:,:) 
    197200      END SELECT 
    198201 
     
    896899 
    897900      !                                !* Check of some namelist values 
    898       IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' ) 
    899       IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' ) 
    900       IF( nn_z0_met < 0 .OR. nn_z0_met > 2 ) CALL ctl_stop( 'bad flag: nn_z0_met is 0, 1 or 2' ) 
    901       IF( nn_stab_func  < 0 .OR. nn_stab_func  > 3 ) CALL ctl_stop( 'bad flag: nn_stab_func is 0, 1, 2 and 3' ) 
    902       IF( nn_clos       < 0 .OR. nn_clos       > 3 ) CALL ctl_stop( 'bad flag: nn_clos is 0, 1, 2 or 3' ) 
     901      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' ) 
     902      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' ) 
     903      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' ) 
     904      IF( nn_z0_met == 3 .AND. .NOT.ln_wave ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_wave=T' ) 
     905      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' ) 
     906      IF( nn_clos       < 0 .OR. nn_clos       > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' ) 
    903907 
    904908      SELECT CASE ( nn_clos )          !* set the parameters for the chosen closure 
  • branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90

    r7359 r7422  
    7070         DO jj = 1, jpjm1 
    7171            DO ji = 1, fs_jpim1 
    72                qbv(ji,jj,jk) = 1.0 * 0.353553 * swh(ji,jj) * tsd2d(ji,jj) *             & 
     72               qbv(ji,jj,jk) = 1.0 * 0.353553 * hsw(ji,jj) * tsd2d(ji,jj) *             & 
    7373            &                  EXP(3.0 * wnum(ji,jj) *                                  &                      
    7474            &                  (-MIN( gdepw_n(ji  ,jj  ,jk), gdepw_n(ji+1,jj  ,jk),     & 
  • branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/step.F90

    r7351 r7422  
    210210                         CALL dyn_adv       ( kstp )  ! advection (vector or flux form) 
    211211                         CALL dyn_vor       ( kstp )  ! vorticity term including Coriolis 
    212       IF( ln_wave .AND. ln_sdw .AND. ln_stcor)           & 
    213                &         CALL dyn_stcor     ( kstp )  ! Stokes-Coriolis forcing 
    214212                         CALL dyn_ldf       ( kstp )  ! lateral mixing 
    215213                         CALL dyn_hpg       ( kstp )  ! horizontal gradient of Hydrostatic pressure 
  • branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

    r7351 r7422  
    4242   USE dynzdf           ! vertical diffusion               (dyn_zdf routine) 
    4343   USE dynspg           ! surface pressure gradient        (dyn_spg routine) 
    44    USE dynstcor         ! simp. form of Stokes-Coriolis 
    4544 
    4645   USE dynnxt           ! time-stepping                    (dyn_nxt routine) 
Note: See TracChangeset for help on using the changeset viewer.