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 11277 for branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90 – NEMO

Ignore:
Timestamp:
2019-07-17T15:29:15+02:00 (5 years ago)
Author:
kingr
Message:

Merged Juan's changes for running AMM15 woth wave coupling.
Corrected minor logic error to allow AMM7-uncoupled to reproduce earlier results.
Few line spacing changes to allow merging with OBS br and trunk rvn 5518.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90

    r8058 r11277  
    44   !! Wave module  
    55   !!====================================================================== 
    6    !! History :  3.3.1  !   2011-09  (Adani M)  Original code: Drag Coefficient  
    7    !!         :  3.4    !   2012-10  (Adani M)                 Stokes Drift  
     6   !! History :  3.3  !   2011-09  (Adani M)  Original code: Drag Coefficient  
     7   !!         :  3.4  !   2012-10  (Adani M)                 Stokes Drift  
     8   !!            3.6  !   2014-09  (Clementi E, Oddo P)New Stokes Drift Computation 
     9   !!             -   !   2016-12  (G. Madec, E. Clementi) update Stoke drift computation  
     10   !!                                                     + add sbc_wave_ini routine 
    811   !!---------------------------------------------------------------------- 
    9    USE iom             ! I/O manager library 
    10    USE in_out_manager  ! I/O manager 
    11    USE lib_mpp         ! distribued memory computing library 
    12    USE fldread        ! read input fields 
    13    USE oce 
    14    USE sbc_oce        ! Surface boundary condition: ocean fields 
    15    USE domvvl 
    16  
    17     
     12 
    1813   !!---------------------------------------------------------------------- 
    19    !!   sbc_wave       : read drag coefficient from wave model in netcdf files  
     14   !!   sbc_stokes    : calculate 3D Stokes-drift velocities 
     15   !!   sbc_wave      : wave data from wave model in netcdf files  
     16   !!   sbc_wave_init : initialisation fo surface waves 
    2017   !!---------------------------------------------------------------------- 
     18   USE oce            ! ocean variables 
     19   USE sbc_oce        ! Surface boundary condition: ocean fields 
     20   USE bdy_oce        ! open boundary condition variables 
     21   USE domvvl         ! domain: variable volume layers 
     22   ! 
     23   USE iom            ! I/O manager library 
     24   USE in_out_manager ! I/O manager 
     25   USE lib_mpp        ! distribued memory computing library 
     26   USE fldread        ! read input fields 
     27   USE wrk_nemo       ! 
     28   USE phycst         ! physical constants  
    2129 
    2230   IMPLICIT NONE 
    2331   PRIVATE 
    2432 
    25    PUBLIC   sbc_wave    ! routine called in sbc_blk_core or sbc_blk_mfs 
     33   PUBLIC   sbc_stokes      ! routine called in sbccpl  
     34   PUBLIC   sbc_stress      ! routine called in sbcmod 
     35   PUBLIC   sbc_wave        ! routine called in sbcmod  
     36   PUBLIC   sbc_wave_init   ! routine called in sbcmod 
    2637    
    27    INTEGER , PARAMETER ::   jpfld  = 3           ! maximum number of files to read for srokes drift 
    28    INTEGER , PARAMETER ::   jp_usd = 1           ! index of stokes drift  (i-component) (m/s)    at T-point 
    29    INTEGER , PARAMETER ::   jp_vsd = 2           ! index of stokes drift  (j-component) (m/s)    at T-point 
    30    INTEGER , PARAMETER ::   jp_wn  = 3           ! index of wave number                 (1/m)    at T-point 
    31    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_cd    ! structure of input fields (file informations, fields read) Drag Coefficient 
    32    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_sd    ! structure of input fields (file informations, fields read) Stokes Drift 
    33    REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: cdn_wave  
    34    REAL(wp),ALLOCATABLE,DIMENSION (:,:)              :: usd2d,vsd2d,uwavenum,vwavenum  
    35    REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:,:)     :: usd3d,vsd3d,wsd3d  
    36  
    37    !! * Substitutions 
    38 #  include "domzgr_substitute.h90" 
     38   ! Variables checking if the wave parameters are coupled (if not, they are read from file) 
     39   LOGICAL, PUBLIC ::   cpl_hsig   = .FALSE. 
     40   LOGICAL, PUBLIC ::   cpl_phioc  = .FALSE. 
     41   LOGICAL, PUBLIC ::   cpl_sdrft  = .FALSE. 
     42   LOGICAL, PUBLIC ::   cpl_wper   = .FALSE. 
     43   LOGICAL, PUBLIC ::   cpl_wfreq  = .FALSE. 
     44   LOGICAL, PUBLIC ::   cpl_wnum   = .FALSE. 
     45   LOGICAL, PUBLIC ::   cpl_tauoc  = .FALSE. 
     46   LOGICAL, PUBLIC ::   cpl_tauw   = .FALSE. 
     47   LOGICAL, PUBLIC ::   cpl_wdrag  = .FALSE. 
     48 
     49   INTEGER ::   nn_sdrift      ! type of parameterization to calculate vertical Stokes drift 
     50   INTEGER, PARAMETER ::   jp_breivik  = 0     ! Breivik 2015: v_z=v_0*[exp(2*k*z)/(1-8*k*z)] 
     51   INTEGER, PARAMETER ::   jp_phillips = 1     ! Phillips:     v_z=v_o*[exp(2*k*z)-beta*sqrt(-2*k*pi*z)*erfc(sqrt(-2*k*z))] 
     52   INTEGER, PARAMETER ::   jp_peakph   = 2     ! Phillips using the peak wave number read from wave model instead of the inverse depth scale 
     53 
     54   INTEGER ::   jpfld    ! number of files to read for stokes drift 
     55   INTEGER ::   jp_usd   ! index of stokes drift  (i-component) (m/s)    at T-point 
     56   INTEGER ::   jp_vsd   ! index of stokes drift  (j-component) (m/s)    at T-point 
     57   INTEGER ::   jp_hsw   ! index of significant wave hight      (m)      at T-point 
     58   INTEGER ::   jp_wmp   ! index of mean wave period            (s)      at T-point 
     59   INTEGER ::   jp_wfr   ! index of wave peak frequency         (s^-1)   at T-point 
     60 
     61   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_cd    ! structure of input fields (file informations, fields read) Drag Coefficient 
     62   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_sd    ! structure of input fields (file informations, fields read) Stokes Drift 
     63   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_wn    ! structure of input fields (file informations, fields read) wave number for Qiao 
     64   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_tauoc ! structure of input fields (file informations, fields read) normalized wave stress into the ocean 
     65   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_tauw ! structure of input fields (file informations, fields read) ocean stress components from wave model 
     66   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_phioc ! structure of input fields (file informations, fields read) wave to ocean energy  
     67   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   cdn_wave            !: 
     68   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   hsw, wmp, wnum      !: 
     69   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   wfreq               !:  
     70   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   rn_crban            !: Craig and Banner constant for surface breaking waves mixing 
     71   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauoc_wave          !: 
     72   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauw_x              !: 
     73   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauw_y              !: 
     74   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tsd2d               !: 
     75   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   div_sd              !: barotropic stokes drift divergence 
     76   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   ut0sd, vt0sd        !: surface Stokes drift velocities at t-point 
     77   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) ::   usd  , vsd  , wsd   !: Stokes drift velocities at u-, v- & w-points, resp. 
     78 
     79#  include "vectopt_loop_substitute.h90" 
    3980   !!---------------------------------------------------------------------- 
    4081   !! NEMO/OPA 4.0 , NEMO Consortium (2011)  
     
    4485CONTAINS 
    4586 
     87   SUBROUTINE sbc_stokes( ) 
     88      !!--------------------------------------------------------------------- 
     89      !!                     ***  ROUTINE sbc_stokes  *** 
     90      !! 
     91      !! ** Purpose :   compute the 3d Stokes Drift according to Breivik et al., 
     92      !!                2014 (DOI: 10.1175/JPO-D-14-0020.1) 
     93      !! 
     94      !! ** Method  : - Calculate Stokes transport speed  
     95      !!              - Calculate horizontal divergence  
     96      !!              - Integrate the horizontal divergenze from the bottom  
     97      !! ** action   
     98      !!--------------------------------------------------------------------- 
     99      INTEGER  ::   jj, ji, jk   ! dummy loop argument 
     100      INTEGER  ::   ik           ! local integer  
     101      REAL(wp) ::  ztransp, zfac, ztemp 
     102      REAL(wp) ::  zdep_u, zdep_v, zkh_u, zkh_v, zda_u, zda_v 
     103      REAL(wp), DIMENSION(:,:)  , POINTER ::   zk_t, zk_u, zk_v, zu0_sd, zv0_sd   ! 2D workspace 
     104      REAL(wp), DIMENSION(:,:,:), POINTER ::   ze3divh                            ! 3D workspace 
     105 
     106      !!--------------------------------------------------------------------- 
     107      ! 
     108 
     109      CALL wrk_alloc( jpi,jpj,jpk,   ze3divh ) 
     110      CALL wrk_alloc( jpi,jpj,       zk_t, zk_u, zk_v, zu0_sd, zv0_sd )  
     111      ! 
     112      ! select parameterization for the calculation of vertical Stokes drift 
     113      ! exp. wave number at t-point 
     114      IF( nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips ) THEN   ! (Eq. (19) in Breivick et al. (2014) ) 
     115         zfac = 2.0_wp * rpi / 16.0_wp 
     116         DO jj = 1, jpj                
     117            DO ji = 1, jpi 
     118               ! Stokes drift velocity estimated from Hs and Tmean 
     119               ztransp = zfac * hsw(ji,jj)*hsw(ji,jj) / MAX( wmp(ji,jj), 0.0000001_wp ) 
     120               ! Stokes surface speed 
     121               tsd2d(ji,jj) = SQRT( ut0sd(ji,jj)*ut0sd(ji,jj) + vt0sd(ji,jj)*vt0sd(ji,jj)) 
     122               ! Wavenumber scale 
     123               zk_t(ji,jj) = ABS( tsd2d(ji,jj) ) / MAX( ABS( 5.97_wp*ztransp ), 0.0000001_wp ) 
     124            END DO 
     125         END DO 
     126         DO jj = 1, jpjm1              ! exp. wave number & Stokes drift velocity at u- & v-points 
     127            DO ji = 1, jpim1 
     128               zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) ) 
     129               zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) ) 
     130               ! 
     131               zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) ) 
     132               zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) ) 
     133            END DO 
     134         END DO 
     135      ELSE IF( nn_sdrift==jp_peakph ) THEN    ! peak wave number calculated from the peak frequency received by the wave model 
     136         DO jj = 1, jpjm1               
     137            DO ji = 1, jpim1 
     138               zk_u(ji,jj) = 0.5_wp * ( wfreq(ji,jj)*wfreq(ji,jj) + wfreq(ji+1,jj)*wfreq(ji+1,jj) ) / grav 
     139               zk_v(ji,jj) = 0.5_wp * ( wfreq(ji,jj)*wfreq(ji,jj) + wfreq(ji,jj+1)*wfreq(ji,jj+1) ) / grav 
     140               ! 
     141               zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) ) 
     142               zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) ) 
     143            END DO 
     144         END DO 
     145      ENDIF 
     146      ! 
     147      !                       !==  horizontal Stokes Drift 3D velocity  ==! 
     148      IF( nn_sdrift==jp_breivik ) THEN 
     149         DO jk = 1, jpkm1 
     150            DO jj = 2, jpjm1 
     151               DO ji = 2, jpim1 
     152                  zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) ) 
     153                  zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) ) 
     154                  !                           
     155                  zkh_u = zk_u(ji,jj) * zdep_u     ! k * depth 
     156                  zkh_v = zk_v(ji,jj) * zdep_v 
     157                  !                                ! Depth attenuation 
     158                  zda_u = EXP( -2.0_wp*zkh_u ) / ( 1.0_wp + 8.0_wp*zkh_u ) 
     159                  zda_v = EXP( -2.0_wp*zkh_v ) / ( 1.0_wp + 8.0_wp*zkh_v ) 
     160                  ! 
     161                  usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk) 
     162                  vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk) 
     163               END DO 
     164            END DO 
     165         END DO 
     166      ELSE IF( nn_sdrift==jp_phillips .OR. nn_sdrift==jp_peakph ) THEN 
     167         DO jk = 1, jpkm1 
     168            DO jj = 2, jpjm1 
     169               DO ji = 2, jpim1 
     170                  zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) ) 
     171                  zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) ) 
     172                  !                           
     173                  zkh_u = zk_u(ji,jj) * zdep_u     ! k * depth 
     174                  zkh_v = zk_v(ji,jj) * zdep_v 
     175                  !                                ! Depth attenuation 
     176                  zda_u = EXP( -2.0_wp*zkh_u ) - SQRT(2.0_wp*rpi*zkh_u) * ERFC(SQRT(2.0_wp*zkh_u)) 
     177                  zda_v = EXP( -2.0_wp*zkh_v ) - SQRT(2.0_wp*rpi*zkh_v) * ERFC(SQRT(2.0_wp*zkh_v)) 
     178                  ! 
     179                  usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk) 
     180                  vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk) 
     181               END DO 
     182            END DO 
     183         END DO 
     184      ENDIF 
     185 
     186      CALL lbc_lnk( usd(:,:,:), 'U', vsd(:,:,:), 'V', -1. ) 
     187      ! 
     188      !                       !==  vertical Stokes Drift 3D velocity  ==! 
     189      ! 
     190      DO jk = 1, jpkm1               ! Horizontal e3*divergence 
     191         DO jj = 2, jpj 
     192            DO ji = fs_2, jpi 
     193               ze3divh(ji,jj,jk) = (  e2u(ji  ,jj) * e3u_n(ji  ,jj,jk) * usd(ji,  jj,jk)    & 
     194                  &                 - e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) * usd(ji-1,jj,jk)    & 
     195                  &                 + e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) * vsd(ji,jj  ,jk)    & 
     196                  &                 - e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vsd(ji,jj-1,jk)  ) * r1_e12t(ji,jj) 
     197            END DO 
     198         END DO 
     199      END DO 
     200      ! 
     201      IF( .NOT. AGRIF_Root() ) THEN 
     202         IF( nbondi ==  1 .OR. nbondi == 2 )   ze3divh(nlci-1,   :  ,:) = 0._wp      ! east 
     203         IF( nbondi == -1 .OR. nbondi == 2 )   ze3divh(  2   ,   :  ,:) = 0._wp      ! west 
     204         IF( nbondj ==  1 .OR. nbondj == 2 )   ze3divh(  :   ,nlcj-1,:) = 0._wp      ! north 
     205         IF( nbondj == -1 .OR. nbondj == 2 )   ze3divh(  :   ,  2   ,:) = 0._wp      ! south 
     206      ENDIF 
     207      ! 
     208      CALL lbc_lnk( ze3divh, 'T', 1. ) 
     209      ! 
     210      IF( .NOT. lk_vvl ) THEN   ;   ik = 1   ! none zero velocity through the sea surface 
     211      ELSE                      ;   ik = 2   ! w=0 at the surface (set one for all in sbc_wave_init) 
     212      ENDIF 
     213      DO jk = jpkm1, ik, -1          ! integrate from the bottom the hor. divergence (NB: at k=jpk w is always zero) 
     214         wsd(:,:,jk) = wsd(:,:,jk+1) - ze3divh(:,:,jk) 
     215      END DO 
     216#if defined key_bdy 
     217      IF( lk_bdy ) THEN 
     218         DO jk = 1, jpkm1 
     219            wsd(:,:,jk) = wsd(:,:,jk) * bdytmask(:,:) 
     220         END DO 
     221      ENDIF 
     222#endif 
     223      !                       !==  Horizontal divergence of barotropic Stokes transport  ==! 
     224      div_sd(:,:) = 0._wp 
     225      DO jk = 1, jpkm1                                 !  
     226        div_sd(:,:) = div_sd(:,:) + ze3divh(:,:,jk) 
     227      END DO 
     228      ! 
     229      CALL iom_put( "ustokes",  usd  ) 
     230      CALL iom_put( "vstokes",  vsd  ) 
     231      CALL iom_put( "wstokes",  wsd  ) 
     232      ! 
     233      CALL wrk_dealloc( jpi,jpj,jpk,   ze3divh ) 
     234      CALL wrk_dealloc( jpi,jpj,       zk_t, zk_u, zk_v, zu0_sd, zv0_sd ) 
     235      ! 
     236   END SUBROUTINE sbc_stokes 
     237 
     238 
     239   SUBROUTINE sbc_stress( ) 
     240      !!--------------------------------------------------------------------- 
     241      !!                     ***  ROUTINE sbc_stress  *** 
     242      !! 
     243      !! ** Purpose :   Updates the ocean momentum modified by waves 
     244      !! 
     245      !! ** Method  : - Calculate u,v components of stress depending on stress 
     246      !!                model  
     247      !!              - Calculate the stress module 
     248      !!              - The wind module is not modified by waves  
     249      !! ** action   
     250      !!--------------------------------------------------------------------- 
     251      INTEGER  ::   jj, ji   ! dummy loop argument 
     252      ! 
     253      IF( ln_tauoc ) THEN 
     254         utau(:,:) = utau(:,:)*tauoc_wave(:,:) 
     255         vtau(:,:) = vtau(:,:)*tauoc_wave(:,:) 
     256         taum(:,:) = taum(:,:)*tauoc_wave(:,:) 
     257      ENDIF 
     258      ! 
     259      IF( ln_tauw ) THEN 
     260         DO jj = 1, jpjm1 
     261            DO ji = 1, jpim1 
     262               ! Stress components at u- & v-points 
     263               utau(ji,jj) = 0.5_wp * ( tauw_x(ji,jj) + tauw_x(ji+1,jj) ) 
     264               vtau(ji,jj) = 0.5_wp * ( tauw_y(ji,jj) + tauw_y(ji,jj+1) ) 
     265               ! 
     266               ! Stress module at t points 
     267               taum(ji,jj) = SQRT( tauw_x(ji,jj)*tauw_x(ji,jj) + tauw_y(ji,jj)*tauw_y(ji,jj) ) 
     268            END DO 
     269         END DO 
     270         CALL lbc_lnk_multi( utau(:,:), 'U', -1. , vtau(:,:), 'V', -1. , taum(:,: ), 'T', -1. ) 
     271      ENDIF 
     272      ! 
     273   END SUBROUTINE sbc_stress 
     274 
     275 
    46276   SUBROUTINE sbc_wave( kt ) 
    47277      !!--------------------------------------------------------------------- 
    48       !!                     ***  ROUTINE sbc_apr  *** 
    49       !! 
    50       !! ** Purpose :   read drag coefficient from wave model  in netcdf files. 
     278      !!                     ***  ROUTINE sbc_wave  *** 
     279      !! 
     280      !! ** Purpose :   read wave parameters from wave model  in netcdf files. 
    51281      !! 
    52282      !! ** Method  : - Read namelist namsbc_wave 
    53283      !!              - Read Cd_n10 fields in netcdf files  
    54284      !!              - Read stokes drift 2d in netcdf files  
    55       !!              - Read wave number      in netcdf files  
    56       !!              - Compute 3d stokes drift using monochromatic 
    57       !! ** action  :    
    58       !!                
    59       !!--------------------------------------------------------------------- 
    60       USE oce,  ONLY : un,vn,hdivn,rotn 
    61       USE divcur 
    62       USE wrk_nemo 
    63 #if defined key_bdy 
    64       USE bdy_oce, ONLY : bdytmask 
    65 #endif 
    66       INTEGER, INTENT( in  ) ::  kt       ! ocean time step 
    67       INTEGER                ::  ierror   ! return error code 
    68       INTEGER                ::  ifpr, jj,ji,jk  
    69       INTEGER                ::   ios     ! Local integer output status for namelist read 
    70       REAL(wp),DIMENSION(:,:,:),POINTER             ::  udummy,vdummy,hdivdummy,rotdummy 
    71       REAL                                          ::  z2dt,z1_2dt 
    72       TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i     ! array of namelist informations on the fields to read 
     285      !!              - Read wave number in netcdf files  
     286      !!              - Compute 3d stokes drift using Breivik et al.,2014 
     287      !!                formulation 
     288      !! ** action   
     289      !!--------------------------------------------------------------------- 
     290      INTEGER, INTENT(in   ) ::   kt   ! ocean time step 
     291      !!--------------------------------------------------------------------- 
     292      ! 
     293      IF( ln_cdgw .AND. .NOT. cpl_wdrag ) THEN     !==  Neutral drag coefficient  ==! 
     294         CALL fld_read( kt, nn_fsbc, sf_cd )             ! read from external forcing 
     295         cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1) 
     296         ! check that the drag coefficient contains proper information even if 
     297         ! the masks do not match - the momentum stress is not masked! 
     298         WHERE( cdn_wave < 0.0 ) cdn_wave = 1.5e-3 
     299         WHERE( cdn_wave > 1.0 ) cdn_wave = 1.5e-3 
     300      ENDIF 
     301 
     302      IF( ln_tauoc .AND. .NOT. cpl_tauoc ) THEN    !==  Wave induced stress  ==! 
     303         CALL fld_read( kt, nn_fsbc, sf_tauoc )          ! read wave norm stress from external forcing 
     304         tauoc_wave(:,:) = sf_tauoc(1)%fnow(:,:,1) 
     305         WHERE( tauoc_wave <   0.0 ) tauoc_wave = 1.0 
     306         WHERE( tauoc_wave > 100.0 ) tauoc_wave = 1.0 
     307      ENDIF 
     308 
     309      IF( ln_tauw .AND. .NOT. cpl_tauw ) THEN      !==  Wave induced stress  ==! 
     310         CALL fld_read( kt, nn_fsbc, sf_tauw )           ! read ocean stress components from external forcing (T grid) 
     311         tauw_x(:,:) = sf_tauw(1)%fnow(:,:,1) 
     312         WHERE( tauw_x < -100.0 ) tauw_x = 0.0 
     313         WHERE( tauw_x >  100.0 ) tauw_x = 0.0 
     314 
     315         tauw_y(:,:) = sf_tauw(2)%fnow(:,:,1) 
     316         WHERE( tauw_y < -100.0 ) tauw_y = 0.0 
     317         WHERE( tauw_y >  100.0 ) tauw_y = 0.0 
     318      ENDIF 
     319 
     320      IF( ln_phioc .AND. .NOT. cpl_phioc ) THEN    !==  Wave to ocean energy  ==! 
     321         CALL fld_read( kt, nn_fsbc, sf_phioc )          ! read wave to ocean energy from external forcing 
     322         rn_crban(:,:) = 29.0 * sf_phioc(1)%fnow(:,:,1)     ! ! Alfa is phioc*sqrt(rau0/zrhoa)  : rau0=water density, zhroa= air density 
     323         WHERE( rn_crban > 1.e8   ) rn_crban = 0.0    !remove first mask mistmatch points, then cap values in case of low friction velocity 
     324         WHERE( rn_crban < 0.0    ) rn_crban = 0.0 
     325         WHERE( rn_crban > 1000.0 ) rn_crban = 1000.0 
     326      ENDIF 
     327 
     328      IF( ln_sdw .OR. ln_rough )  THEN             !==  Computation of the 3d Stokes Drift  ==!  
     329         ! 
     330         IF( jpfld > 0 ) THEN                            ! Read from file only if the field is not coupled 
     331            CALL fld_read( kt, nn_fsbc, sf_sd )          ! read wave parameters from external forcing 
     332            IF( jp_hsw > 0 ) THEN 
     333               hsw  (:,:) = sf_sd(jp_hsw)%fnow(:,:,1)   ! significant wave height 
     334               WHERE( hsw > 100.0 ) hsw = 0.0 
     335               WHERE( hsw <   0.0 ) hsw = 0.0 
     336            ENDIF 
     337            IF( jp_wmp > 0 ) THEN 
     338               wmp  (:,:) = sf_sd(jp_wmp)%fnow(:,:,1)   ! wave mean period 
     339               WHERE( wmp > 100.0 ) wmp = 0.0 
     340               WHERE( wmp <   0.0 ) wmp = 0.0 
     341            ENDIF 
     342            IF( jp_wfr > 0 ) THEN 
     343               wfreq(:,:) = sf_sd(jp_wfr)%fnow(:,:,1)   ! Peak wave frequency  
     344               WHERE( wfreq <    0.0 ) wfreq = 0.0  
     345               WHERE( wfreq >  100.0 ) wfreq = 0.0 
     346            ENDIF 
     347            IF( jp_usd > 0 ) THEN 
     348               ut0sd(:,:) = sf_sd(jp_usd)%fnow(:,:,1)   ! 2D zonal Stokes Drift at T point 
     349               WHERE( ut0sd < -100.0 ) ut0sd = 1.0 
     350               WHERE( ut0sd >  100.0 ) ut0sd = 1.0 
     351            ENDIF 
     352            IF( jp_vsd > 0 ) THEN 
     353               vt0sd(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)   ! 2D meridional Stokes Drift at T point 
     354               WHERE( vt0sd < -100.0 ) vt0sd = 1.0 
     355               WHERE( vt0sd >  100.0 ) vt0sd = 1.0 
     356            ENDIF 
     357         ENDIF 
     358      ENDIF 
     359      ! 
     360      IF( ln_sdw ) THEN 
     361         ! Read also wave number if needed, so that it is available in coupling routines 
     362         IF( ln_zdfqiao .AND. .NOT.cpl_wnum ) THEN 
     363            CALL fld_read( kt, nn_fsbc, sf_wn )          ! read wave parameters from external forcing 
     364            wnum(:,:) = sf_wn(1)%fnow(:,:,1) 
     365         ENDIF 
     366            
     367         !                                         !==  Computation of the 3d Stokes Drift  ==!  
     368         ! 
     369         IF( ((nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) .AND. & 
     370                          jp_hsw>0 .AND. jp_wmp>0 .AND. jp_usd>0 .AND. jp_vsd>0) .OR. & 
     371             (nn_sdrift==jp_peakph .AND. jp_wfr>0 .AND. jp_usd>0 .AND. jp_vsd>0) ) & 
     372            CALL sbc_stokes()            ! Calculate only if required fields are read 
     373         !                               ! In coupled wave model-NEMO case the call is done after coupling 
     374         ! 
     375      ENDIF 
     376      ! 
     377   END SUBROUTINE sbc_wave 
     378 
     379 
     380   SUBROUTINE sbc_wave_init 
     381      !!--------------------------------------------------------------------- 
     382      !!                     ***  ROUTINE sbc_wave_init  *** 
     383      !! 
     384      !! ** Purpose :   read wave parameters from wave model  in netcdf files. 
     385      !! 
     386      !! ** Method  : - Read namelist namsbc_wave 
     387      !!              - Read Cd_n10 fields in netcdf files  
     388      !!              - Read stokes drift 2d in netcdf files  
     389      !!              - Read wave number in netcdf files  
     390      !!              - Compute 3d stokes drift using Breivik et al.,2014 
     391      !!                formulation 
     392      !! ** action   
     393      !!--------------------------------------------------------------------- 
     394      INTEGER ::   ierror, ios   ! local integer 
     395      INTEGER ::   ifpr 
     396      !! 
    73397      CHARACTER(len=100)     ::  cn_dir                          ! Root directory for location of drag coefficient files 
    74       TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd, sn_wn   ! informations about the fields to be read 
    75       !!--------------------------------------------------------------------- 
    76       NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_wn 
    77       !!--------------------------------------------------------------------- 
    78  
    79       !!---------------------------------------------------------------------- 
    80       ! 
    81       ! 
    82       !                                         ! -------------------- ! 
    83       IF( kt == nit000 ) THEN                   ! First call kt=nit000 ! 
    84          !                                      ! -------------------- ! 
    85          REWIND( numnam_ref )              ! Namelist namsbc_wave in reference namelist : File for drag coeff. from wave model 
    86          READ  ( numnam_ref, namsbc_wave, IOSTAT = ios, ERR = 901) 
    87 901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in reference namelist', lwp ) 
    88  
    89          REWIND( numnam_cfg )              ! Namelist namsbc_wave in configuration namelist : File for drag coeff. from wave model 
    90          READ  ( numnam_cfg, namsbc_wave, IOSTAT = ios, ERR = 902 ) 
    91 902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in configuration namelist', lwp ) 
    92          IF(lwm) WRITE ( numond, namsbc_wave ) 
    93          ! 
    94  
    95          IF ( ln_cdgw ) THEN 
     398      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i, slf_j     ! array of namelist informations on the fields to read 
     399      TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd, sn_phioc, & 
     400                             &   sn_hsw, sn_wmp, sn_wfr, sn_wnum , & 
     401                             &   sn_tauoc, sn_tauwx, sn_tauwy      ! information about the fields to be read 
     402      ! 
     403      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wfr, sn_wnum, sn_phioc,          & 
     404                             sn_tauoc, sn_tauwx, sn_tauwy,                                                       & 
     405                             ln_cdgw, ln_sdw, ln_stcor, ln_phioc, ln_tauoc, ln_tauw, ln_zdfqiao, ln_rough,       & 
     406                             nn_sdrift, nn_wmix 
     407      !!--------------------------------------------------------------------- 
     408      ! 
     409      REWIND( numnam_ref )              ! Namelist namsbc_wave in reference namelist : File for drag coeff. from wave model 
     410      READ  ( numnam_ref, namsbc_wave, IOSTAT = ios, ERR = 901) 
     411901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in reference namelist', lwp ) 
     412          
     413      REWIND( numnam_cfg )              ! Namelist namsbc_wave in configuration namelist : File for drag coeff. from wave model 
     414      READ  ( numnam_cfg, namsbc_wave, IOSTAT = ios, ERR = 902 ) 
     415902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in configuration namelist', lwp ) 
     416      IF(lwm) WRITE ( numond, namsbc_wave ) 
     417      ! 
     418      IF(lwp) THEN               !* Parameter print 
     419         WRITE(numout,*) 
     420         WRITE(numout,*) 'sbc_wave_init: wave physics' 
     421         WRITE(numout,*) '~~~~~~~~' 
     422         WRITE(numout,*) '   Namelist namsbc_wave : set wave physics parameters' 
     423         WRITE(numout,*) '      Stokes drift corr. to vert. velocity    ln_sdw      = ', ln_sdw  
     424         WRITE(numout,*) '        vertical parametrization              nn_sdrift   = ', nn_sdrift  
     425         WRITE(numout,*) '      Stokes coriolis term                    ln_stcor    = ', ln_stcor  
     426         WRITE(numout,*) '      wave modified ocean stress              ln_tauoc    = ', ln_tauoc  
     427         WRITE(numout,*) '      wave modified ocean stress components   ln_tauw     = ', ln_tauw  
     428         WRITE(numout,*) '      wave to ocean energy                    ln_phioc    = ', ln_phioc 
     429         WRITE(numout,*) '        vertical mixing parametrization       nn_wmix     = ', nn_wmix  
     430         WRITE(numout,*) '      neutral drag coefficient                ln_cdgw     = ', ln_cdgw 
     431         WRITE(numout,*) '      wave roughness length modification      ln_rough    = ', ln_rough  
     432         WRITE(numout,*) '      Qiao vertical mixing formulation        ln_zdfqiao  = ', ln_zdfqiao 
     433      ENDIF 
     434 
     435      IF ( ln_wave ) THEN 
     436         ! Activated wave physics but no wave physics components activated  
     437         IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_tauw .OR. ln_stcor .OR. ln_phioc & 
     438                                                                    .OR. ln_rough .OR. ln_zdfqiao) )   THEN   
     439             CALL ctl_warn( 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauoc=F, ln_tauw=F, ln_stcor=F ', & 
     440                                                      'ln_phioc=F, ln_rough=F, ln_zdfqiao=F' )  
     441         ELSE 
     442         IF (ln_stcor .AND. .NOT. ln_sdw) &  
     443             CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T') 
     444         IF ( ln_cdgw .AND. .NOT.(nn_drag==jp_ukmo .OR. nn_drag==jp_std .OR. nn_drag==jp_const .OR. nn_drag==jp_mcore) ) &  
     445             CALL ctl_stop( 'The chosen nn_drag for momentum calculation must be 0, 1, 2, or 3') 
     446         IF ( ln_cdgw .AND. ln_blk_core .AND. nn_drag==0 ) & 
     447             CALL ctl_stop( 'The chosen nn_drag for momentum calculation in core forcing must be 1, 2, or 3') 
     448         IF ( ln_cdgw .AND. ln_flx .AND. nn_drag==3 ) & 
     449             CALL ctl_stop( 'The chosen nn_drag for momentum calculation in direct forcing must be 0, 1, or 2') 
     450         IF( ln_phioc .AND. .NOT.(nn_wmix==jp_craigbanner .OR. nn_wmix==jp_janssen) ) &  
     451            CALL ctl_stop( 'The chosen nn_wmix for wave vertical mixing must be 0, or 1' ) 
     452         IF( ln_sdw .AND. .NOT.(nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips .OR. nn_sdrift==jp_peakph) ) &  
     453            CALL ctl_stop( 'The chosen nn_sdrift for Stokes drift vertical velocity must be 0, 1, or 2' ) 
     454         IF( ln_zdfqiao .AND. .NOT.ln_sdw ) &  
     455            CALL ctl_stop( 'Qiao vertical mixing can not be used without Stokes drift (ln_sdw)' ) 
     456         IF( ln_tauoc .AND. ln_tauw ) &  
     457            CALL ctl_stop( 'More than one method for modifying the ocean stress has been selected ', & 
     458                                     '(ln_tauoc=.true. and ln_tauw=.true.)' ) 
     459         IF( ln_tauoc ) & 
     460             CALL ctl_warn( 'You are subtracting the wave stress to the ocean (ln_tauoc=.true.)' ) 
     461         IF( ln_tauw ) & 
     462             CALL ctl_warn( 'The wave modified ocean stress components are used (ln_tauw=.true.) ', & 
     463                                  'This will override any other specification of the ocean stress' )  
     464         ENDIF 
     465      ELSE 
     466         IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_tauw .OR. ln_stcor .OR. ln_phioc .OR. ln_rough .OR. ln_zdfqiao ) &   
     467            &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    &  
     468            &                  'with drag coefficient (ln_cdgw =T) '  ,                        &  
     469            &                  'or Stokes Drift (ln_sdw=T) ' ,                                 &  
     470            &                  'or Stokes-Coriolis term (ln_stcor=T)',                         & 
     471            &                  'or ocean stress modification due to waves (ln_tauoc=T) ',      &  
     472            &                  'or ocean stress components from waves (ln_tauw=T) ',          &  
     473            &                  'or wave to ocean energy modification (ln_phioc=T) ',           &  
     474            &                  'or wave surface roughness (ln_rough=T) ',                      &  
     475            &                  'or Qiao vertical mixing formulation (ln_zdfqiao=T) ' ) 
     476      ENDIF  
     477      ! 
     478      IF( ln_cdgw ) THEN 
     479         IF( .NOT. cpl_wdrag ) THEN 
    96480            ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg 
    97             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
     481            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_cd structure' ) 
    98482            ! 
    99483                                   ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1)   ) 
    100484            IF( sn_cdg%ln_tint )   ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) ) 
    101             CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 
    102             ALLOCATE( cdn_wave(jpi,jpj) ) 
    103             cdn_wave(:,:) = 0.0 
    104         ENDIF 
    105          IF ( ln_sdw ) THEN 
    106             slf_i(jp_usd) = sn_usd ; slf_i(jp_vsd) = sn_vsd; slf_i(jp_wn) = sn_wn 
    107             ALLOCATE( sf_sd(3), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg 
    108             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
     485            CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' ) 
     486         ENDIF 
     487         ALLOCATE( cdn_wave(jpi,jpj) ) 
     488      ENDIF 
     489 
     490      IF( ln_tauoc ) THEN 
     491         IF( .NOT. cpl_tauoc ) THEN 
     492            ALLOCATE( sf_tauoc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_tauoc 
     493            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_tauoc structure' ) 
    109494            ! 
    110             DO ifpr= 1, jpfld 
    111                ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) ) 
    112                IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) ) 
    113             END DO 
    114             CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 
    115             ALLOCATE( usd2d(jpi,jpj),vsd2d(jpi,jpj),uwavenum(jpi,jpj),vwavenum(jpi,jpj) ) 
    116             ALLOCATE( usd3d(jpi,jpj,jpk),vsd3d(jpi,jpj,jpk),wsd3d(jpi,jpj,jpk) ) 
    117             usd2d(:,:) = 0.0 ;  vsd2d(:,:) = 0.0 ; uwavenum(:,:) = 0.0 ; vwavenum(:,:) = 0.0 
    118             usd3d(:,:,:) = 0.0 ;vsd3d(:,:,:) = 0.0 ; wsd3d(:,:,:) = 0.0 
    119          ENDIF 
    120       ENDIF 
     495                                    ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1)   ) 
     496            IF( sn_tauoc%ln_tint )  ALLOCATE( sf_tauoc(1)%fdta(jpi,jpj,1,2) ) 
     497            CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' ) 
     498         ENDIF 
     499         ALLOCATE( tauoc_wave(jpi,jpj) ) 
     500      ENDIF 
     501 
     502      IF( ln_tauw ) THEN 
     503         IF( .NOT. cpl_tauw ) THEN 
     504            ALLOCATE( sf_tauw(2), STAT=ierror )           !* allocate and fill sf_wave with sn_tauwx/y 
     505            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_tauw structure' ) 
     506            ! 
     507            ALLOCATE( slf_j(2) ) 
     508            slf_j(1) = sn_tauwx 
     509            slf_j(2) = sn_tauwy 
     510                                    ALLOCATE( sf_tauw(1)%fnow(jpi,jpj,1)   ) 
     511                                    ALLOCATE( sf_tauw(2)%fnow(jpi,jpj,1)   ) 
     512            IF( slf_j(1)%ln_tint )  ALLOCATE( sf_tauw(1)%fdta(jpi,jpj,1,2) ) 
     513            IF( slf_j(2)%ln_tint )  ALLOCATE( sf_tauw(2)%fdta(jpi,jpj,1,2) ) 
     514            CALL fld_fill( sf_tauw, (/ slf_j /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' ) 
     515         ENDIF 
     516         ALLOCATE( tauw_x(jpi,jpj) ) 
     517         ALLOCATE( tauw_y(jpi,jpj) ) 
     518      ENDIF 
     519 
     520      IF( ln_phioc ) THEN 
     521         IF( .NOT. cpl_phioc ) THEN 
     522            ALLOCATE( sf_phioc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_phioc 
     523            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_phioc structure' ) 
     524            ! 
     525                                    ALLOCATE( sf_phioc(1)%fnow(jpi,jpj,1)   ) 
     526            IF( sn_phioc%ln_tint )  ALLOCATE( sf_phioc(1)%fdta(jpi,jpj,1,2) ) 
     527            CALL fld_fill( sf_phioc, (/ sn_phioc /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' ) 
     528         ENDIF 
     529         ALLOCATE( rn_crban(jpi,jpj) ) 
     530      ENDIF 
     531 
     532      ! Find out how many fields have to be read from file if not coupled 
     533      jpfld=0 
     534      jp_usd=0   ;   jp_vsd=0   ;   jp_hsw=0   ;   jp_wmp=0   ;   jp_wfr=0 
     535      IF( ln_sdw ) THEN 
     536         IF( .NOT. cpl_sdrft ) THEN 
     537            jpfld  = jpfld + 1 
     538            jp_usd = jpfld 
     539            jpfld  = jpfld + 1 
     540            jp_vsd = jpfld 
     541         ENDIF 
     542         IF( .NOT. cpl_hsig .AND. (nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) ) THEN 
     543            jpfld  = jpfld + 1 
     544            jp_hsw = jpfld 
     545         ENDIF 
     546         IF( .NOT. cpl_wper .AND. (nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) ) THEN 
     547            jpfld  = jpfld + 1 
     548            jp_wmp = jpfld 
     549         ENDIF 
     550         IF( .NOT. cpl_wfreq .AND. nn_sdrift==jp_peakph ) THEN 
     551            jpfld  = jpfld + 1 
     552            jp_wfr = jpfld 
     553         ENDIF 
     554      ENDIF 
     555 
     556      IF( ln_rough .AND. .NOT. cpl_hsig .AND. jp_hsw==0 ) THEN 
     557         jpfld  = jpfld + 1 
     558         jp_hsw = jpfld 
     559      ENDIF 
     560 
     561      ! Read from file only the non-coupled fields  
     562      IF( jpfld > 0 ) THEN 
     563         ALLOCATE( slf_i(jpfld) ) 
     564         IF( jp_usd > 0 )   slf_i(jp_usd) = sn_usd 
     565         IF( jp_vsd > 0 )   slf_i(jp_vsd) = sn_vsd 
     566         IF( jp_hsw > 0 )   slf_i(jp_hsw) = sn_hsw 
     567         IF( jp_wmp > 0 )   slf_i(jp_wmp) = sn_wmp 
     568         IF( jp_wfr > 0 )   slf_i(jp_wfr) = sn_wfr 
     569         ALLOCATE( sf_sd(jpfld), STAT=ierror )   !* allocate and fill sf_sd with stokes drift 
     570         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_sd structure' ) 
    121571         ! 
     572         DO ifpr= 1, jpfld 
     573            ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) ) 
     574            IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) ) 
     575         END DO 
    122576         ! 
    123       IF ( ln_cdgw ) THEN 
    124          CALL fld_read( kt, nn_fsbc, sf_cd )      !* read drag coefficient from external forcing 
    125          cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1) 
    126       ENDIF 
    127       IF ( ln_sdw )  THEN 
    128           CALL fld_read( kt, nn_fsbc, sf_sd )      !* read drag coefficient from external forcing 
    129  
    130          ! Interpolate wavenumber, stokes drift into the grid_V and grid_V 
    131          !------------------------------------------------- 
    132  
    133          DO jj = 1, jpjm1 
    134             DO ji = 1, jpim1 
    135                uwavenum(ji,jj)=0.5 * ( 2. - umask(ji,jj,1) ) * ( sf_sd(3)%fnow(ji,jj,1) * tmask(ji,jj,1) & 
    136                &                                + sf_sd(3)%fnow(ji+1,jj,1) * tmask(ji+1,jj,1) ) 
    137  
    138                vwavenum(ji,jj)=0.5 * ( 2. - vmask(ji,jj,1) ) * ( sf_sd(3)%fnow(ji,jj,1) * tmask(ji,jj,1) & 
    139                &                                + sf_sd(3)%fnow(ji,jj+1,1) * tmask(ji,jj+1,1) ) 
    140  
    141                usd2d(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( sf_sd(1)%fnow(ji,jj,1) * tmask(ji,jj,1) & 
    142                &                                + sf_sd(1)%fnow(ji+1,jj,1) * tmask(ji+1,jj,1) ) 
    143  
    144                vsd2d(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( sf_sd(2)%fnow(ji,jj,1) * tmask(ji,jj,1) & 
    145                &                                + sf_sd(2)%fnow(ji,jj+1,1) * tmask(ji,jj+1,1) ) 
    146             END DO 
    147          END DO 
    148  
    149           !Computation of the 3d Stokes Drift 
    150           DO jk = 1, jpk 
    151              DO jj = 1, jpj-1 
    152                 DO ji = 1, jpi-1 
    153                    usd3d(ji,jj,jk) = usd2d(ji,jj)*exp(2.0*uwavenum(ji,jj)*(-MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj  ,jk)))) 
    154                    vsd3d(ji,jj,jk) = vsd2d(ji,jj)*exp(2.0*vwavenum(ji,jj)*(-MIN( gdept_0(ji,jj,jk) , gdept_0(ji  ,jj+1,jk)))) 
    155                 END DO 
    156              END DO 
    157              usd3d(jpi,:,jk) = usd2d(jpi,:)*exp( 2.0*uwavenum(jpi,:)*(-gdept_0(jpi,:,jk)) ) 
    158              vsd3d(:,jpj,jk) = vsd2d(:,jpj)*exp( 2.0*vwavenum(:,jpj)*(-gdept_0(:,jpj,jk)) ) 
    159           END DO 
    160  
    161           CALL wrk_alloc( jpi,jpj,jpk,udummy,vdummy,hdivdummy,rotdummy) 
    162            
    163           udummy(:,:,:)=un(:,:,:) 
    164           vdummy(:,:,:)=vn(:,:,:) 
    165           hdivdummy(:,:,:)=hdivn(:,:,:) 
    166           rotdummy(:,:,:)=rotn(:,:,:) 
    167           un(:,:,:)=usd3d(:,:,:) 
    168           vn(:,:,:)=vsd3d(:,:,:) 
    169           CALL div_cur(kt) 
    170       !                                           !------------------------------! 
    171       !                                           !     Now Vertical Velocity    ! 
    172       !                                           !------------------------------! 
    173           z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog) 
    174  
    175           z1_2dt = 1.e0 / z2dt 
    176           DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence 
    177              ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise 
    178              wsd3d(:,:,jk) = wsd3d(:,:,jk+1) -   fse3t_n(:,:,jk) * hdivn(:,:,jk)        & 
    179                 &                      - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    & 
    180                 &                         * tmask(:,:,jk) * z1_2dt 
    181 #if defined key_bdy 
    182              wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:) 
    183 #endif 
    184           END DO 
    185           hdivn(:,:,:)=hdivdummy(:,:,:) 
    186           rotn(:,:,:)=rotdummy(:,:,:) 
    187           vn(:,:,:)=vdummy(:,:,:) 
    188           un(:,:,:)=udummy(:,:,:) 
    189           CALL wrk_dealloc( jpi,jpj,jpk,udummy,vdummy,hdivdummy,rotdummy) 
    190       ENDIF 
    191    END SUBROUTINE sbc_wave 
    192        
     577         CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' ) 
     578      ENDIF 
     579 
     580      IF( ln_sdw ) THEN 
     581         ALLOCATE( usd  (jpi,jpj,jpk), vsd  (jpi,jpj,jpk), wsd(jpi,jpj,jpk) ) 
     582         ALLOCATE( wmp  (jpi,jpj)  ) 
     583         ALLOCATE( wfreq (jpi,jpj) ) 
     584         ALLOCATE( ut0sd(jpi,jpj)    , vt0sd(jpi,jpj)     ) 
     585         ALLOCATE( div_sd(jpi,jpj) ) 
     586         ALLOCATE( tsd2d (jpi,jpj) ) 
     587         usd(:,:,:) = 0._wp 
     588         vsd(:,:,:) = 0._wp 
     589         wsd(:,:,:) = 0._wp 
     590         ! Wave number needed only if ln_zdfqiao=T 
     591         IF( ln_zdfqiao .AND. .NOT.cpl_wnum ) THEN 
     592            ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fill sf_wave with sn_wnum 
     593            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable toallocate sf_wn structure' ) 
     594                                   ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1)   ) 
     595            IF( sn_wnum%ln_tint )  ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) ) 
     596            CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'read wave input', 'namsbc_wave' ) 
     597         ENDIF 
     598         ALLOCATE( wnum(jpi,jpj) ) 
     599      ENDIF 
     600 
     601      IF( ln_sdw .OR. ln_rough ) THEN 
     602         ALLOCATE( hsw  (jpi,jpj) ) 
     603      ENDIF 
     604      ! 
     605   END SUBROUTINE sbc_wave_init 
     606 
    193607   !!====================================================================== 
    194608END MODULE sbcwave 
Note: See TracChangeset for help on using the changeset viewer.