Changeset 9115


Ignore:
Timestamp:
2017-12-18T15:06:42+01:00 (3 years ago)
Author:
acc
Message:

Branch 2017/dev_merge_2017. Tidier version of sbcwave.F90 with associated changes in modules. References and schemes for Stokes drift parameterisations have been updated. Choices now are the original Breivik 2014 scheme and the latest Li et al 2017 scheme (which is based on the Breivik et al 2016 Phillips spectrum but with a depth averaged profile). This has been compiled but not yet tested. Also removed trailing spaces in lbc_lnk_multi_generic.h90

Location:
branches/2017/dev_merge_2017
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_merge_2017/DOC/TexFiles/Bibliography/Biblio.bib

    r9052 r9115  
    545545URL = {http://www.geosci-model-dev.net/8/1285/2015/}, 
    546546DOI = {10.5194/gmd-8-1285-2015} 
     547} 
     548 
     549@ARTICLE{Breivik_al_OM2016, 
     550   AUTHOR = {{\O}yvind Breivik and Jean-Raymond Bidlot and Peter A.E.M. Janssen}, 
     551   YEAR = {2016}, 
     552   TITLE = "{A Stokes drift approximation based on the Phillips spectrum}", 
     553   JOURNAL = {OM}, 
     554   VOLUME = {100}, 
     555   DOI = {10.1016/j.ocemod.2016.01.005}, 
     556   PAGES = {49--56, arXiv:1601.08092} 
    547557} 
    548558 
     
    21222132  author = {S Levitus }, 
    21232133  pages = {173 pp} 
     2134} 
     2135 
     2136@ARTICLE{Li_al_OM2017}, 
     2137   AUTHOR = {Q Li and B Fox-Kemper and {\O} Breivik and A Webb}, 
     2138   YEAR = {2017}, 
     2139   TITLE = {Statistical Models of Global Langmuir Mixing},  
     2140   JOURNAL = {OM}, 
     2141   VOLUME = {113}, 
     2142   ISSUE = {May}, 
     2143   PAGES = {95--114}, 
     2144   DOI = {10.1016/j.ocemod.2017.03.016} 
    21242145} 
    21252146 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LBC/lbc_lnk_multi_generic.h90

    r9093 r9115  
    1 #if defined DIM_2d  
    2 #   define ARRAY_TYPE(i,j,k,l)   REAL(wp), DIMENSION(i,j)  
    3 #   define PTR_TYPE              TYPE(PTR_2D)  
    4 #   define PTR_ptab              pt2d  
     1#if defined DIM_2d 
     2#   define ARRAY_TYPE(i,j,k,l)   REAL(wp), DIMENSION(i,j) 
     3#   define PTR_TYPE              TYPE(PTR_2D) 
     4#   define PTR_ptab              pt2d 
    55#endif 
    66#if defined DIM_3d 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r9033 r9115  
    8383   !!           Stokes drift parametrization definition  
    8484   !!---------------------------------------------------------------------- 
    85    INTEGER , PUBLIC, PARAMETER ::   jp_breivik  = 0     ! Breivik 2015: v_z=v_0*[exp(2*k*z)/(1-8*k*z)] 
    86    INTEGER , PUBLIC, PARAMETER ::   jp_phillips = 1     ! Phillips: v_z=v_o*[exp(2*k*z)-beta*sqrt(-2*k*pi*z)*erfc(sqrt(-2*k*z))] 
    87    INTEGER , PUBLIC, PARAMETER ::   jp_peakfr   = 2     ! Phillips using the peak wave number read from wave model instead of the inverse depth scale 
     85   INTEGER , PUBLIC, PARAMETER ::   jp_breivik_2014 = 0     !: Breivik  2014: v_z=v_0*[exp(2*k*z)/(1-8*k*z)] 
     86   INTEGER , PUBLIC, PARAMETER ::   jp_li_2017      = 1     !: Li et al 2017: Stokes drift based on Phillips spectrum (Breivik 2016)  
     87                                                            !  with depth averaged profile 
     88   INTEGER , PUBLIC, PARAMETER ::   jp_peakfr       = 2     !: Li et al 2017: using the peak wave number read from wave model instead  
     89                                                            !  of the inverse depth scale 
     90   LOGICAL , PUBLIC            ::   ll_st_bv2014  = .FALSE. !  logical indicator, .true. if Breivik 2014 parameterisation is active. 
     91   LOGICAL , PUBLIC            ::   ll_st_li2017  = .FALSE. !  logical indicator, .true. if Li 2017 parameterisation is active. 
     92   LOGICAL , PUBLIC            ::   ll_st_bv_li   = .FALSE. !  logical indicator, .true. if either Breivik or Li parameterisation is active. 
     93   LOGICAL , PUBLIC            ::   ll_st_peakfr  = .FALSE. !  logical indicator, .true. if using Li 2017 with peak wave number 
    8894 
    8995   !!---------------------------------------------------------------------- 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r9033 r9115  
    163163      ! 
    164164      IF( ln_sdw ) THEN 
    165          IF( .NOT.(nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips .OR. nn_sdrift==jp_peakfr) ) & 
     165         IF( .NOT.(nn_sdrift==jp_breivik_2014 .OR. nn_sdrift==jp_li_2017 .OR. nn_sdrift==jp_peakfr) ) & 
    166166            CALL ctl_stop( 'The chosen nn_sdrift for Stokes drift vertical velocity must be 0, 1, or 2' ) 
    167167      ENDIF 
     168      ll_st_bv2014  = ( nn_sdrift==jp_breivik_2014 ) 
     169      ll_st_li2017  = ( nn_sdrift==jp_li_2017 ) 
     170      ll_st_bv_li   = ( ll_st_bv2014 .OR. ll_st_li2017 ) 
     171      ll_st_peakfr  = ( nn_sdrift==jp_peakfr ) 
    168172      IF( ln_tauwoc .AND. ln_tauw ) & 
    169173         CALL ctl_stop( 'More than one method for modifying the ocean stress has been selected ', & 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90

    r9098 r9115  
    4545   LOGICAL, PUBLIC ::   cpl_wfreq  = .FALSE. 
    4646   LOGICAL, PUBLIC ::   cpl_wnum   = .FALSE. 
    47    LOGICAL, PUBLIC ::   cpl_tauwoc  = .FALSE. 
     47   LOGICAL, PUBLIC ::   cpl_tauwoc = .FALSE. 
    4848   LOGICAL, PUBLIC ::   cpl_tauw   = .FALSE. 
    4949   LOGICAL, PUBLIC ::   cpl_wdrag  = .FALSE. 
     
    5959   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_sd      ! structure of input fields (file informations, fields read) Stokes Drift 
    6060   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_wn      ! structure of input fields (file informations, fields read) wave number for Qiao 
    61    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tauwoc   ! structure of input fields (file informations, fields read) normalized wave stress into the ocean 
     61   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tauwoc  ! structure of input fields (file informations, fields read) normalized wave stress into the ocean 
    6262   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tauw    ! structure of input fields (file informations, fields read) ocean stress components from wave model 
    6363 
     
    100100      REAL(wp) ::  zstokes_psi_u_bot, zstokes_psi_v_bot 
    101101      REAL(wp) ::  zdep_u, zdep_v, zkh_u, zkh_v 
    102       REAL(wp), DIMENSION(:,:)  , POINTER ::   zk_t, zk_u, zk_v, zu0_sd, zv0_sd   ! 2D workspace 
    103       REAL(wp), DIMENSION(:,:)  , POINTER ::   zstokes_psi_u_top, zstokes_psi_v_top   ! 2D workspace 
    104       REAL(wp), DIMENSION(:,:,:), POINTER ::   ze3divh                            ! 3D workspace 
    105       !!--------------------------------------------------------------------- 
    106       ! 
    107       CALL wrk_alloc( jpi,jpj,jpk,   ze3divh ) 
    108       CALL wrk_alloc( jpi,jpj,       zk_t, zk_u, zk_v, zu0_sd, zv0_sd ) 
    109       CALL wrk_alloc( jpi,jpj,       zstokes_psi_u_top, zstokes_psi_v_top) 
     102      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   zk_t, zk_u, zk_v, zu0_sd, zv0_sd     ! 2D workspace 
     103      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   zstokes_psi_u_top, zstokes_psi_v_top ! 2D workspace 
     104      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ze3divh                              ! 3D workspace 
     105      !!--------------------------------------------------------------------- 
     106      ! 
     107      ALLOCATE( ze3divh(jpi,jpj,jpk) ) 
     108      ALLOCATE( zk_t(jpi,jpj), zk_u(jpi,jpj), zk_v(jpi,jpj), zu0_sd(jpi,jpj), zv0_sd(jpi,jpj) ) 
    110109      ! 
    111110      ! select parameterization for the calculation of vertical Stokes drift 
    112111      ! exp. wave number at t-point 
    113       IF( nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips ) THEN   ! (Eq. (19) in Breivick et al. (2014) ) 
     112      IF( ll_st_bv_li ) THEN   ! (Eq. (19) in Breivik et al. (2014) ) 
    114113         zfac = 2.0_wp * rpi / 16.0_wp 
    115114         DO jj = 1, jpj 
     
    132131            END DO 
    133132         END DO 
    134       ELSE IF( nn_sdrift==jp_peakfr ) THEN    ! peak wave number calculated from the peak frequency received by the wave model 
     133      ELSE IF( ll_st_peakfr ) THEN    ! peak wave number calculated from the peak frequency received by the wave model 
    135134         DO jj = 1, jpjm1 
    136135            DO ji = 1, jpim1 
     
    145144      ! 
    146145      !                       !==  horizontal Stokes Drift 3D velocity  ==! 
    147       IF( nn_sdrift==jp_breivik ) THEN 
     146      IF( ll_st_bv2014 ) THEN 
    148147         DO jk = 1, jpkm1 
    149148            DO jj = 2, jpjm1 
     
    163162            END DO 
    164163         END DO 
    165       ELSE IF( nn_sdrift==jp_phillips .OR. nn_sdrift==jp_peakfr ) THEN 
     164      ELSE IF( ll_st_li2017 .OR. ll_st_peakfr ) THEN 
     165         ALLOCATE( zstokes_psi_u_top(jpi,jpj), zstokes_psi_v_top(jpi,jpj) ) 
     166         DO jj = 1, jpjm1              ! exp. wave number & Stokes drift velocity at u- & v-points 
     167            DO ji = 1, jpim1 
     168               zstokes_psi_u_top(ji,jj) = 0._wp 
     169               zstokes_psi_v_top(ji,jj) = 0._wp 
     170            END DO 
     171         END DO 
    166172         DO jk = 1, jpkm1 
    167173            DO jj = 2, jpjm1 
    168174               DO ji = 2, jpim1 
    169                   zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) ) 
    170                   zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) ) 
    171                   !                           
    172                   zkh_u = zk_u(ji,jj) * zdep_u     ! k * depth 
    173                   zkh_v = zk_v(ji,jj) * zdep_v 
    174                   !                                ! Depth attenuation 
    175                   zda_u = EXP( -2.0_wp*zkh_u ) - SQRT(2.0_wp*rpi*zkh_u) * ERFC(SQRT(2.0_wp*zkh_u)) 
    176                   zda_v = EXP( -2.0_wp*zkh_v ) - SQRT(2.0_wp*rpi*zkh_v) * ERFC(SQRT(2.0_wp*zkh_v)) 
     175                  zbot_u = 0.5_wp * ( gdepw_n(ji,jj,jk+1) + gdepw_n(ji+1,jj,jk+1) ) 
     176                  zbot_v = 0.5_wp * ( gdepw_n(ji,jj,jk+1) + gdepw_n(ji,jj+1,jk+1) ) 
     177                  zkb_u  = 2.0_wp * zk_u(ji,jj) * zbot_u                             ! 2k * bottom depth 
     178                  zkb_v  = 2.0_wp * zk_v(ji,jj) * zbot_v                             ! 2k * bottom depth 
     179                  ! 
     180                  zke3_u = MAX(1.e-8_wp, 2.0_wp * zk_u(ji,jj) * e3u_n(ji,jj,jk))     ! 2k * thickness 
     181                  zke3_v = MAX(1.e-8_wp, 2.0_wp * zk_v(ji,jj) * e3v_n(ji,jj,jk))     ! 2k * thickness 
     182 
     183                  ! Depth attenuation .... do u component first.. 
     184                  zdepth      = zkb_u 
     185                  zsqrt_depth = SQRT(zdepth) 
     186                  zexp_depth  = EXP(-zdepth) 
     187                  zstokes_psi_u_bot = 1.0_wp - zexp_depth  & 
     188                       &              - z_two_thirds * ( zsqrtpi*zsqrt_depth*zdepth*ERFC(zsqrt_depth) & 
     189                       &              + 1.0_wp - (1.0_wp + zdepth)*zexp_depth ) 
     190                  zda_u                    = ( zstokes_psi_u_bot - zstokes_psi_u_top(ji,jj) ) / zke3_u 
     191                  zstokes_psi_u_top(ji,jj) =   zstokes_psi_u_bot 
     192 
     193                  !         ... and then v component 
     194                  zdepth      =zkb_v 
     195                  zsqrt_depth = SQRT(zdepth) 
     196                  zexp_depth  = EXP(-zdepth) 
     197                  zstokes_psi_v_bot = 1.0_wp - zexp_depth  & 
     198                       &              - z_two_thirds * ( zsqrtpi*zsqrt_depth*zdepth*ERFC(zsqrt_depth) & 
     199                       &              + 1.0_wp - (1.0_wp + zdepth)*zexp_depth ) 
     200                  zda_v                    = ( zstokes_psi_v_bot - zstokes_psi_v_top(ji,jj) ) / zke3_v 
     201                  zstokes_psi_v_top(ji,jj) =   zstokes_psi_v_bot 
    177202                  ! 
    178203                  usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk) 
     
    181206            END DO 
    182207         END DO 
     208         DEALLOCATE( zstokes_psi_u_top, zstokes_psi_v_top ) 
    183209      ENDIF 
    184210 
     
    232258      CALL iom_put( "wstokes",  wsd  ) 
    233259      ! 
    234       CALL wrk_dealloc( jpi,jpj,jpk,   ze3divh ) 
    235       CALL wrk_dealloc( jpi,jpj,       zk_t, zk_u, zk_v, zu0_sd, zv0_sd ) 
    236       CALL wrk_dealloc( jpi,jpj,       zstokes_psi_u_top, zstokes_psi_v_top) 
     260      DEALLOCATE( ze3divh ) 
     261      DEALLOCATE( zk_t, zk_u, zk_v, zu0_sd, zv0_sd ) 
    237262      ! 
    238263   END SUBROUTINE sbc_stokes 
     
    298323      ENDIF 
    299324 
    300       IF( ln_tauwoc .AND. .NOT. cpl_tauwoc ) THEN    !==  Wave induced stress  ==! 
    301          CALL fld_read( kt, nn_fsbc, sf_tauwoc )          ! read wave norm stress from external forcing 
     325      IF( ln_tauwoc .AND. .NOT. cpl_tauwoc ) THEN  !==  Wave induced stress  ==! 
     326         CALL fld_read( kt, nn_fsbc, sf_tauwoc )         ! read wave norm stress from external forcing 
    302327         tauoc_wave(:,:) = sf_tauwoc(1)%fnow(:,:,1) 
    303328      ENDIF 
     
    326351         ENDIF 
    327352            
    328          !                                         !==  Computation of the 3d Stokes Drift  ==!  
     353         ! Calculate only if required fields have been read 
     354         ! In coupled wave model-NEMO case the call is done after coupling 
    329355         ! 
    330          IF( ((nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) .AND. & 
    331                           jp_hsw>0 .AND. jp_wmp>0 .AND. jp_usd>0 .AND. jp_vsd>0) .OR. & 
    332              (nn_sdrift==jp_peakfr .AND. jp_wfr>0 .AND. jp_usd>0 .AND. jp_vsd>0) ) & 
    333             CALL sbc_stokes()            ! Calculate only if required fields are read 
    334          !                               ! In coupled wave model-NEMO case the call is done after coupling 
     356         IF( ( ll_st_bv_li   .AND. jp_hsw>0 .AND. jp_wmp>0 .AND. jp_usd>0 .AND. jp_vsd>0 ) .OR. & 
     357           & ( ll_st_peakfr  .AND. jp_wfr>0 .AND. jp_usd>0 .AND. jp_vsd>0                ) ) CALL sbc_stokes() 
    335358         ! 
    336359      ENDIF 
     
    356379      INTEGER ::   ifpr 
    357380      !! 
    358       CHARACTER(len=100)     ::  cn_dir                          ! Root directory for location of drag coefficient files 
     381      CHARACTER(len=100)     ::  cn_dir                            ! Root directory for location of drag coefficient files 
    359382      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i, slf_j     ! array of namelist informations on the fields to read 
    360383      TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd,  & 
    361384                             &   sn_hsw, sn_wmp, sn_wfr, sn_wnum, & 
    362                              &   sn_tauwoc, sn_tauwx, sn_tauwy      ! informations about the fields to be read 
     385                             &   sn_tauwoc, sn_tauwx, sn_tauwy     ! informations about the fields to be read 
    363386      ! 
    364387      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wfr, & 
     
    377400      IF( ln_cdgw ) THEN 
    378401         IF( .NOT. cpl_wdrag ) THEN 
    379             ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg 
     402            ALLOCATE( sf_cd(1), STAT=ierror )               !* allocate and fill sf_wave with sn_cdg 
    380403            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' ) 
    381404            ! 
     
    392415            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' ) 
    393416            ! 
    394                                     ALLOCATE( sf_tauwoc(1)%fnow(jpi,jpj,1)   ) 
     417                                     ALLOCATE( sf_tauwoc(1)%fnow(jpi,jpj,1)   ) 
    395418            IF( sn_tauwoc%ln_tint )  ALLOCATE( sf_tauwoc(1)%fdta(jpi,jpj,1,2) ) 
    396419            CALL fld_fill( sf_tauwoc, (/ sn_tauwoc /), cn_dir, 'sbc_wave_init', 'Wave module', 'namsbc_wave' ) 
     
    428451            jp_vsd = jpfld 
    429452         ENDIF 
    430          IF( .NOT. cpl_hsig .AND. (nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) ) THEN 
     453         IF( .NOT. cpl_hsig  .AND. ll_st_bv_li ) THEN 
    431454            jpfld  = jpfld + 1 
    432455            jp_hsw = jpfld 
    433456         ENDIF 
    434          IF( .NOT. cpl_wper .AND. (nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) ) THEN 
     457         IF( .NOT. cpl_wper  .AND. ll_st_bv_li ) THEN 
    435458            jpfld  = jpfld + 1 
    436459            jp_wmp = jpfld 
    437460         ENDIF 
    438          IF( .NOT. cpl_wfreq .AND. nn_sdrift==jp_peakfr ) THEN 
     461         IF( .NOT. cpl_wfreq .AND. ll_st_peakfr ) THEN 
    439462            jpfld  = jpfld + 1 
    440463            jp_wfr = jpfld 
Note: See TracChangeset for help on using the changeset viewer.