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 8750 for branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90 – NEMO

Ignore:
Timestamp:
2017-11-20T12:45:31+01:00 (6 years ago)
Author:
jcastill
Message:

First set of changes for ticket #1980
Addition of the Phillips parametrization for vertical Stokes drift

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90

    r8749 r8750  
    4242   LOGICAL, PUBLIC ::   cpl_sdrfty = .FALSE. 
    4343   LOGICAL, PUBLIC ::   cpl_wper   = .FALSE. 
     44   LOGICAL, PUBLIC ::   cpl_wfreq  = .FALSE. 
    4445   LOGICAL, PUBLIC ::   cpl_wnum   = .FALSE. 
    4546   LOGICAL, PUBLIC ::   cpl_wstrf  = .FALSE. 
     
    5152   INTEGER ::   jp_hsw   ! index of significant wave hight      (m)      at T-point 
    5253   INTEGER ::   jp_wmp   ! index of mean wave period            (s)      at T-point 
     54   INTEGER ::   jp_wfr   ! index of wave peak frequency         (1/s)    at T-point 
    5355 
    5456   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_cd      ! structure of input fields (file informations, fields read) Drag Coefficient 
     
    5860   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   cdn_wave            !: 
    5961   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   hsw, wmp, wnum      !:  
     62   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   wfreq               !:  
    6063   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauoc_wave          !:   
    6164   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tsd2d               !:  
     
    9699      CALL wrk_alloc( jpi,jpj,       zk_t, zk_u, zk_v, zu0_sd, zv0_sd ) 
    97100      ! 
    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 
     101      ! select parameterization for the calculation of vertical Stokes drift 
     102      ! exp. wave number at t-point 
     103      IF( nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips ) THEN   ! (Eq. (19) in Breivick et al. (2014) ) 
     104         zfac = 2.0_wp * rpi / 16.0_wp 
     105         DO jj = 1, jpj 
     106            DO ji = 1, jpi 
    102107               ! Stokes drift velocity estimated from Hs and Tmean 
    103                ztransp = zfac * hsw(ji,jj)*hsw(ji,jj) / MAX( wmp(ji,jj) , 0.0000001_wp ) 
     108               ztransp = zfac * hsw(ji,jj)*hsw(ji,jj) / MAX( wmp(ji,jj), 0.0000001_wp ) 
    104109               ! Stokes surface speed 
    105                zsp0 = SQRT( ut0sd(ji,jj)*ut0sd(ji,jj) + vt0sd(ji,jj)*vt0sd(ji,jj) ) 
    106                tsd2d(ji,jj) = zsp0 
     110               tsd2d(ji,jj) = SQRT( ut0sd(ji,jj)*ut0sd(ji,jj) + vt0sd(ji,jj)*vt0sd(ji,jj)) 
    107111               ! Wavenumber scale 
    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 
     112               zk_t(ji,jj) = ABS( tsd2d(ji,jj) ) / MAX( ABS( 5.97_wp*ztransp ), 0.0000001_wp ) 
     113            END DO 
     114         END DO 
     115         DO jj = 1, jpjm1              ! exp. wave number & Stokes drift velocity at u- & v-points 
     116            DO ji = 1, jpim1 
     117               zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) ) 
     118               zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) ) 
     119               ! 
     120               zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) ) 
     121               zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) ) 
     122            END DO 
     123         END DO 
     124      ELSE IF( nn_sdrift==jp_peakfr ) THEN    ! peak wave number calculated from the peak frequency received by the wave model 
     125         DO jj = 1, jpjm1 
     126            DO ji = 1, jpim1 
     127               zk_u(ji,jj) = 0.5_wp * ( wfreq(ji,jj)*wfreq(ji,jj) + wfreq(ji+1,jj)*wfreq(ji+1,jj) ) / grav 
     128               zk_v(ji,jj) = 0.5_wp * ( wfreq(ji,jj)*wfreq(ji,jj) + wfreq(ji,jj+1)*wfreq(ji,jj+1) ) / grav 
     129               ! 
     130               zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) ) 
     131               zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) ) 
     132            END DO 
     133         END DO 
     134      ENDIF 
    120135      ! 
    121136      !                       !==  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 ) 
    133                ! 
    134                usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk) 
    135                vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk) 
    136             END DO 
    137          END DO 
    138       END DO    
     137      IF( nn_sdrift==jp_breivik ) THEN 
     138         DO jk = 1, jpkm1 
     139            DO jj = 2, jpjm1 
     140               DO ji = 2, jpim1 
     141                  zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) ) 
     142                  zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) ) 
     143                  !                           
     144                  zkh_u = zk_u(ji,jj) * zdep_u     ! k * depth 
     145                  zkh_v = zk_v(ji,jj) * zdep_v 
     146                  !                                ! Depth attenuation 
     147                  zda_u = EXP( -2.0_wp*zkh_u ) / ( 1.0_wp + 8.0_wp*zkh_u ) 
     148                  zda_v = EXP( -2.0_wp*zkh_v ) / ( 1.0_wp + 8.0_wp*zkh_v ) 
     149                  ! 
     150                  usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk) 
     151                  vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk) 
     152               END DO 
     153            END DO 
     154         END DO 
     155      ELSE IF( nn_sdrift==jp_phillips .OR. nn_sdrift==jp_peakfr ) THEN 
     156         DO jk = 1, jpkm1 
     157            DO jj = 2, jpjm1 
     158               DO ji = 2, jpim1 
     159                  zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) ) 
     160                  zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) ) 
     161                  !                           
     162                  zkh_u = zk_u(ji,jj) * zdep_u     ! k * depth 
     163                  zkh_v = zk_v(ji,jj) * zdep_v 
     164                  !                                ! Depth attenuation 
     165                  zda_u = EXP( -2.0_wp*zkh_u ) - SQRT(2.0_wp*rpi*zkh_u) * ERFC(SQRT(2.0_wp*zkh_u)) 
     166                  zda_v = EXP( -2.0_wp*zkh_v ) - SQRT(2.0_wp*rpi*zkh_v) * ERFC(SQRT(2.0_wp*zkh_v)) 
     167                  ! 
     168                  usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk) 
     169                  vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk) 
     170               END DO 
     171            END DO 
     172         END DO 
     173      ENDIF 
     174 
    139175      CALL lbc_lnk( usd(:,:,:), 'U', vsd(:,:,:), 'V', -1. ) 
    140176      ! 
     
    222258            IF( jp_hsw > 0 )   hsw  (:,:) = sf_sd(jp_hsw)%fnow(:,:,1)   ! significant wave height 
    223259            IF( jp_wmp > 0 )   wmp  (:,:) = sf_sd(jp_wmp)%fnow(:,:,1)   ! wave mean period 
     260            IF( jp_wfr > 0 )   wfreq(:,:) = sf_sd(jp_wfr)%fnow(:,:,1)   ! Peak wave frequency 
    224261            IF( jp_usd > 0 )   ut0sd(:,:) = sf_sd(jp_usd)%fnow(:,:,1)   ! 2D zonal Stokes Drift at T point 
    225262            IF( jp_vsd > 0 )   vt0sd(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)   ! 2D meridional Stokes Drift at T point 
     
    234271         !                                         !==  Computation of the 3d Stokes Drift  ==!  
    235272         ! 
    236          IF( jpfld == 4 )   CALL sbc_stokes()            ! Calculate only if required fields are read 
    237          !                                               ! In coupled wave model-NEMO case the call is done after coupling 
     273         IF( ((nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) .AND. & 
     274                          jp_hsw>0 .AND. jp_wmp>0 .AND. jp_usd>0 .AND. jp_vsd>0) .OR. & 
     275             (nn_sdrift==jp_peakfr .AND. jp_wfr>0 .AND. jp_usd>0 .AND. jp_vsd>0) ) & 
     276            CALL sbc_stokes()            ! Calculate only if required fields are read 
     277         !                               ! In coupled wave model-NEMO case the call is done after coupling 
    238278         ! 
    239279      ENDIF 
     
    262302      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i     ! array of namelist informations on the fields to read 
    263303      TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd,  & 
    264                              &   sn_hsw, sn_wmp, sn_wnum, sn_tauoc      ! informations about the fields to be read 
    265       ! 
    266       NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wnum, sn_tauoc 
     304                             &   sn_hsw, sn_wmp, sn_wfr, sn_wnum, sn_tauoc      ! informations about the fields to be read 
     305      ! 
     306      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wfr, sn_wnum, sn_tauoc 
    267307      !!--------------------------------------------------------------------- 
    268308      ! 
     
    302342      IF( ln_sdw ) THEN   ! Find out how many fields have to be read from file if not coupled 
    303343         jpfld=0 
    304          jp_usd=0   ;   jp_vsd=0   ;   jp_hsw=0   ;   jp_wmp=0 
     344         jp_usd=0   ;   jp_vsd=0   ;   jp_hsw=0   ;   jp_wmp=0   ;   jp_wfr=0 
    305345         IF( .NOT. cpl_sdrftx ) THEN 
    306346            jpfld  = jpfld + 1 
     
    311351            jp_vsd = jpfld 
    312352         ENDIF 
    313          IF( .NOT. cpl_hsig ) THEN 
     353         IF( .NOT. cpl_hsig .AND. (nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) ) THEN 
    314354            jpfld  = jpfld + 1 
    315355            jp_hsw = jpfld 
    316356         ENDIF 
    317          IF( .NOT. cpl_wper ) THEN 
     357         IF( .NOT. cpl_wper .AND. (nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) ) THEN 
    318358            jpfld  = jpfld + 1 
    319359            jp_wmp = jpfld 
     360         ENDIF 
     361         IF( .NOT. cpl_wfreq .AND. nn_sdrift==jp_peakfr ) THEN 
     362            jpfld  = jpfld + 1 
     363            jp_wfr = jpfld 
    320364         ENDIF 
    321365 
     
    327371            IF( jp_hsw > 0 )   slf_i(jp_hsw) = sn_hsw 
    328372            IF( jp_wmp > 0 )   slf_i(jp_wmp) = sn_wmp 
     373            IF( jp_wfr > 0 )   slf_i(jp_wfr) = sn_wfr 
     374 
    329375            ALLOCATE( sf_sd(jpfld), STAT=ierror )   !* allocate and fill sf_sd with stokes drift 
    330376            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' ) 
     
    339385         ALLOCATE( usd  (jpi,jpj,jpk), vsd  (jpi,jpj,jpk), wsd(jpi,jpj,jpk) ) 
    340386         ALLOCATE( hsw  (jpi,jpj)    , wmp  (jpi,jpj)     ) 
     387         ALLOCATE( wfreq(jpi,jpj) ) 
    341388         ALLOCATE( ut0sd(jpi,jpj)    , vt0sd(jpi,jpj)     ) 
    342389         ALLOCATE( div_sd(jpi,jpj) ) 
Note: See TracChangeset for help on using the changeset viewer.