Changeset 4824


Ignore:
Timestamp:
2014-10-24T18:13:17+02:00 (6 years ago)
Author:
poddo
Message:

#1404: surface wave module

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4822_INGV_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90

    r4624 r4824  
    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.6    !2014-09  (Clementi E, Oddo P)New Stokes Drift Computation  
     7   !!         :  3.3.1  !2011-09  (Adani M)  Original code: Drag Coefficient  
     8   !!         :  3.4    !2012-10  (Adani M)                 Stokes Drift  
    89   !!---------------------------------------------------------------------- 
    910   USE iom             ! I/O manager library 
    1011   USE in_out_manager  ! I/O manager 
    1112   USE lib_mpp         ! distribued memory computing library 
    12    USE fldread        ! read input fields 
     13   USE fldread         ! read input fields 
    1314   USE oce 
    14    USE sbc_oce        ! Surface boundary condition: ocean fields 
     15   USE sbc_oce         ! Surface boundary condition: ocean fields 
    1516   USE domvvl 
    16  
     17   USE phycst 
    1718    
    1819   !!---------------------------------------------------------------------- 
     
    2526   PUBLIC   sbc_wave    ! routine called in sbc_blk_core or sbc_blk_mfs 
    2627    
    27    INTEGER , PARAMETER ::   jpfld  = 3           ! maximum number of files to read for srokes drift 
     28   INTEGER , PARAMETER ::   jpfld  = 4           ! number of files to read for stokes drift 
    2829   INTEGER , PARAMETER ::   jp_usd = 1           ! index of stokes drift  (i-component) (m/s)    at T-point 
    2930   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 
     31   INTEGER , PARAMETER ::   jp_swh = 3           ! index of significant wave hight      (m)      at T-point 
     32   INTEGER , PARAMETER ::   jp_wmp = 4           ! index of mean wave period            (s)      at T-point 
     33! 
     34   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_cd    ! structure of input fields (file informations, fields read) Drag Coefficient 
     35   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_sd    ! structure of input fields (file informations, fields read) Stokes Drift 
     36   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_wn    ! structure of input fields (file informations, fields read) wave number for Qiao 
    3337   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: cdn_wave  
    34    REAL(wp),ALLOCATABLE,DIMENSION (:,:)              :: usd2d,vsd2d,uwavenum,vwavenum  
     38   REAL(wp),ALLOCATABLE,DIMENSION (:,:)              :: usd2d,vsd2d 
     39   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: swh,wmp,wnum 
     40   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: usd2dt,vsd2dt,tsd2d 
    3541   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:,:)     :: usd3d,vsd3d,wsd3d  
    3642 
     
    5460      !!              - Read stokes drift 2d in netcdf files  
    5561      !!              - Read wave number      in netcdf files  
    56       !!              - Compute 3d stokes drift using monochromatic 
     62      !!              - Compute 3d stokes drift using Breivik et al.,2014 
     63      !!                formulation 
    5764      !! ** action  :    
    5865      !!                
     
    6168      USE divcur 
    6269      USE wrk_nemo 
     70      USE zdf_oce,  ONLY : ln_zdfqiao 
    6371#if defined key_bdy 
    6472      USE bdy_oce, ONLY : bdytmask 
     
    6977      INTEGER                ::   ios     ! Local integer output status for namelist read 
    7078      REAL(wp),DIMENSION(:,:,:),POINTER             ::  udummy,vdummy,hdivdummy,rotdummy 
    71       REAL                                          ::  z2dt,z1_2dt 
     79      REAL(wp)                                      ::  z2dt,z1_2dt 
     80      REAL(wp)                                      ::  ztransp,zsp0, zk, zus,zvs 
     81      REAL(wp), DIMENSION(jpi,jpj)                  ::  zfac  
    7282      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i     ! array of namelist informations on the fields to read 
    7383      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 
     84      TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd,  & 
     85                             &   sn_swh, sn_wmp, sn_wnum         ! informations about the fields to be read 
     86      !!--------------------------------------------------------------------- 
     87      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_swh, sn_wmp, sn_wnum 
    7788      !!--------------------------------------------------------------------- 
    7889 
     
    103114            cdn_wave(:,:) = 0.0 
    104115        ENDIF 
     116 
    105117         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 
     118            slf_i(jp_usd) = sn_usd ; slf_i(jp_vsd) = sn_vsd; 
     119            slf_i(jp_swh) = sn_swh ; slf_i(jp_wmp) = sn_wmp; 
     120            ALLOCATE( sf_sd(jpfld), STAT=ierror )           !* allocate and fill sf_sd with stokes drift 
    108121            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
    109122            ! 
     
    112125               IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) ) 
    113126            END DO 
     127! 
    114128            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) ) 
     129            ALLOCATE( usd2d(jpi,jpj),vsd2d(jpi,jpj),usd2dt(jpi,jpj),vsd2dt(jpi,jpj)) 
    116130            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 
     131            ALLOCATE( swh(jpi,jpj), wmp(jpi,jpj) ) 
     132            usd2d(:,:)   = 0.0 ;  vsd2d(:,:)   = 0.0 ; 
     133            usd2dt(:,:)  = 0.0 ;  vsd2dt(:,:)  = 0.0 ; 
     134            usd3d(:,:,:) = 0.0 ;  vsd3d(:,:,:) = 0.0 ; 
     135            wsd3d(:,:,:) = 0.0 ; 
     136            swh(:,:)     = 0.0 ;  wmp(:,:)     = 0.0 ; 
     137            IF ( ln_zdfqiao ) THEN 
     138               ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fillsf_wave with sn_wnum 
     139               IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable toallocate sf_wave structure' ) 
     140                                      ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1)   ) 
     141               IF( sn_wnum%ln_tint )  ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) ) 
     142               CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' ) 
     143               ALLOCATE( wnum(jpi,jpj),tsd2d(jpi,jpj) ) 
     144               wnum(:,:) = 0.0; tsd2d(:,:) = 0.0 
     145            ENDIF 
    119146         ENDIF 
    120147      ENDIF 
    121          ! 
    122148         ! 
    123149      IF ( ln_cdgw ) THEN 
     
    126152      ENDIF 
    127153      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 
     154          CALL fld_read( kt, nn_fsbc, sf_sd )      !* read wave parameters from external forcing 
     155          swh(:,:)    = sf_sd(jp_swh)%fnow(:,:,1) 
     156          wmp(:,:)    = sf_sd(jp_wmp)%fnow(:,:,1) 
     157          usd2dt(:,:) = sf_sd(jp_usd)%fnow(:,:,1)  
     158          vsd2dt(:,:) = sf_sd(jp_vsd)%fnow(:,:,1) 
    131159         !------------------------------------------------- 
    132  
     160         ! Interpolate stokes drift into the grid_V and grid_V 
     161         !------------------------------------------------- 
    133162         DO jj = 1, jpjm1 
    134163            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) ) 
     164               usd2d(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( usd2dt(ji,jj) * tmask(ji,jj,1) & 
     165               &                                             +  usd2dt(ji+1,jj) * tmask(ji+1,jj,1) ) 
     166 
     167               vsd2d(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( vsd2dt(ji,jj) * tmask(ji,jj,1) & 
     168               &                                              + vsd2dt(ji,jj+1) * tmask(ji,jj+1,1) ) 
    146169            END DO 
    147170         END DO 
    148  
    149           !Computation of the 3d Stokes Drift 
     171          !Computation of the 3d Stokes Drift according to Breivik et al.,2014 
     172          !(DOI: 10.1175/JPO-D-14-0020.1) 
    150173          DO jk = 1, jpk 
     174             DO jj = 1, jpj 
     175                DO ji = 1, jpi 
     176 
     177               ! On T grid 
     178               ! Stokes transport speed estimated from Hs and Tmean 
     179               ztransp = 2.0_wp*rpi*swh(ji,jj)**2.0_wp/(16.0_wp*MAX(wmp(ji,jj),0.0000001_wp)) 
     180 
     181               ! Stokes surface speed 
     182               zsp0 = SQRT( sf_sd(jp_usd)%fnow(ji,jj,1)**2 +  sf_sd(jp_vsd)%fnow(ji,jj,1)**2) 
     183 
     184               ! Wavenumber scale 
     185               zk = ABS(zsp0)/MAX(ABS(5.97_wp*ztransp),0.0000001_wp) 
     186 
     187               ! Depth attenuation 
     188               zfac(ji,jj) = EXP(-2.0_wp*zk*fsdept(ji,jj,jk))/(1.0_wp+8.0_wp*zk*fsdept(ji,jj,jk)) 
     189 
     190                END DO 
     191             END DO 
     192! 
    151193             DO jj = 1, jpj-1 
    152194                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)))) 
     195                 ! Into the U and V Grid  
     196                 zus = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zfac(ji,jj) * tmask(ji,jj,1) & 
     197                 &                                + zfac(ji+1,jj) * tmask(ji+1,jj,1) ) 
     198 
     199                 zvs = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zfac(ji,jj) * tmask(ji,jj,1) & 
     200                 &                                + zfac(ji,jj+1) * tmask(ji,jj+1,1) ) 
     201 
     202                 usd3d(ji,jj,jk) = usd2d(ji,jj)*zus 
     203                 vsd3d(ji,jj,jk) = vsd2d(ji,jj)*zvs 
    155204                END DO 
    156205             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)) ) 
     206 
    159207          END DO 
    160208 
    161209          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(:,:,:) 
     210      !-------------------------------------------------------------------           
     211      ! Store horiz. velocity divergence and Rot in TMP array  
     212      ! ------------------------------------------------------------------- 
     213 
     214          udummy(:,:,:)     = un(:,:,:) 
     215          vdummy(:,:,:)     = vn(:,:,:) 
     216          hdivdummy(:,:,:)  = hdivn(:,:,:) 
     217          rotdummy(:,:,:)   = rotn(:,:,:) 
     218          un(:,:,:)         = usd3d(:,:,:) 
     219          vn(:,:,:)         = vsd3d(:,:,:) 
     220 
     221      ! Compute divergence using 3d stokes drift to calculate vertical SD 
     222      ! velocity 
     223      ! ------------------------------------------------------------------- 
     224 
    169225          CALL div_cur(kt) 
     226 
     227      ! ------------------------------------------------------------------- 
    170228      !                                           !------------------------------! 
    171229      !                                           !     Now Vertical Velocity    ! 
     
    177235             ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise 
    178236             wsd3d(:,:,jk) = wsd3d(:,:,jk+1) -   fse3t_n(:,:,jk) * hdivn(:,:,jk)        & 
    179                 &                      - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    & 
    180                 &                         * tmask(:,:,jk) * z1_2dt 
     237                &                            - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    & 
     238                &                                * tmask(:,:,jk) * z1_2dt 
    181239#if defined key_bdy 
    182240             wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:) 
    183241#endif 
    184242          END DO 
    185           hdivn(:,:,:)=hdivdummy(:,:,:) 
    186           rotn(:,:,:)=rotdummy(:,:,:) 
    187           vn(:,:,:)=vdummy(:,:,:) 
    188           un(:,:,:)=udummy(:,:,:) 
     243      ! Back to state variables   
     244      ! ------------------------------------------------------------------- 
     245          hdivn(:,:,:)  = hdivdummy(:,:,:) 
     246          rotn(:,:,:)   = rotdummy(:,:,:) 
     247          vn(:,:,:)     = vdummy(:,:,:) 
     248          un(:,:,:)     = udummy(:,:,:) 
     249 
    189250          CALL wrk_dealloc( jpi,jpj,jpk,udummy,vdummy,hdivdummy,rotdummy) 
     251! 
     252        IF ( ln_zdfqiao )  THEN 
     253          wnum(:,:) = sf_wn(1)%fnow(:,:,1) 
     254         ! Calculate the module of the stokes drift on T grid 
     255         !------------------------------------------------- 
     256         DO jj = 1, jpj 
     257            DO ji = 1, jpi 
     258                tsd2d(ji,jj) = ((sf_sd(jp_usd)%fnow(ji,jj,1) * tmask(ji,jj,1))**2.0  +     & 
     259                &               (sf_sd(jp_vsd)%fnow(ji,jj,1) * tmask(ji,jj,1))**2.0)**0.5 
     260            END DO 
     261         END DO 
     262        ENDIF 
    190263      ENDIF 
     264 
     265 
    191266   END SUBROUTINE sbc_wave 
    192267       
Note: See TracChangeset for help on using the changeset viewer.