- Timestamp:
- 2017-12-13T18:08:50+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
r9019 r9023 33 33 34 34 PUBLIC sbc_stokes ! routine called in sbccpl 35 PUBLIC sbc_wstress ! routine called in sbcmod 35 36 PUBLIC sbc_wave ! routine called in sbcmod 36 37 PUBLIC sbc_wave_init ! routine called in sbcmod … … 42 43 LOGICAL, PUBLIC :: cpl_sdrfty = .FALSE. 43 44 LOGICAL, PUBLIC :: cpl_wper = .FALSE. 45 LOGICAL, PUBLIC :: cpl_wfreq = .FALSE. 44 46 LOGICAL, PUBLIC :: cpl_wnum = .FALSE. 45 LOGICAL, PUBLIC :: cpl_wstrf = .FALSE. 47 LOGICAL, PUBLIC :: cpl_tauoc = .FALSE. 48 LOGICAL, PUBLIC :: cpl_tauw = .FALSE. 46 49 LOGICAL, PUBLIC :: cpl_wdrag = .FALSE. 47 50 … … 51 54 INTEGER :: jp_hsw ! index of significant wave hight (m) at T-point 52 55 INTEGER :: jp_wmp ! index of mean wave period (s) at T-point 56 INTEGER :: jp_wfr ! index of wave peak frequency (1/s) at T-point 53 57 54 58 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_cd ! structure of input fields (file informations, fields read) Drag Coefficient … … 56 60 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_wn ! structure of input fields (file informations, fields read) wave number for Qiao 57 61 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_tauoc ! structure of input fields (file informations, fields read) normalized wave stress into the ocean 62 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_tauw ! structure of input fields (file informations, fields read) ocean stress components from wave model 63 58 64 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: cdn_wave !: 59 65 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: hsw, wmp, wnum !: 66 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: wfreq !: 60 67 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: tauoc_wave !: 68 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: tauw_x, tauw_y !: 61 69 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: tsd2d !: 62 70 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: div_sd !: barotropic stokes drift divergence … … 100 108 CALL wrk_alloc( jpi,jpj, zstokes_psi_u_top, zstokes_psi_v_top) 101 109 ! 102 ! 103 zsqrtpi = SQRT(rpi)104 z_two_thirds = 2.0_wp / 3.0_wp105 zfac =2.0_wp * rpi / 16.0_wp106 DO jj = 1, jpj ! exp. wave number at t-point (Eq. (19) in Breivick et al. (2014) )107 DO ji = 1, jpi110 ! select parameterization for the calculation of vertical Stokes drift 111 ! exp. wave number at t-point 112 IF( nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips ) THEN ! (Eq. (19) in Breivick et al. (2014) ) 113 zfac = 2.0_wp * rpi / 16.0_wp 114 DO jj = 1, jpj 115 DO ji = 1, jpi 108 116 ! Stokes drift velocity estimated from Hs and Tmean 109 ztransp = zfac * hsw(ji,jj)*hsw(ji,jj) / MAX( wmp(ji,jj) 117 ztransp = zfac * hsw(ji,jj)*hsw(ji,jj) / MAX( wmp(ji,jj), 0.0000001_wp ) 110 118 ! Stokes surface speed 111 zsp0 = SQRT( ut0sd(ji,jj)*ut0sd(ji,jj) + vt0sd(ji,jj)*vt0sd(ji,jj) ) 112 tsd2d(ji,jj) = zsp0 119 tsd2d(ji,jj) = SQRT( ut0sd(ji,jj)*ut0sd(ji,jj) + vt0sd(ji,jj)*vt0sd(ji,jj)) 113 120 ! Wavenumber scale 114 zk_t(ji,jj) = (1.0_wp-2.0_wp/3.0_wp)*zsp0/MAX(2.0_wp*ztransp,0.0000001_wp) 115 END DO 116 END DO 117 ! 118 DO jj = 1, jpjm1 ! exp. wave number & Stokes drift velocity at u- & v-points 119 DO ji = 1, jpim1 120 zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) ) 121 zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) ) 122 ! 123 zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) ) 124 zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) ) 125 END DO 126 END DO 121 zk_t(ji,jj) = ABS( tsd2d(ji,jj) ) / MAX( ABS( 5.97_wp*ztransp ), 0.0000001_wp ) 122 END DO 123 END DO 124 DO jj = 1, jpjm1 ! exp. wave number & Stokes drift velocity at u- & v-points 125 DO ji = 1, jpim1 126 zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) ) 127 zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) ) 128 ! 129 zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) ) 130 zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) ) 131 END DO 132 END DO 133 ELSE IF( nn_sdrift==jp_peakfr ) THEN ! peak wave number calculated from the peak frequency received by the wave model 134 DO jj = 1, jpjm1 135 DO ji = 1, jpim1 136 zk_u(ji,jj) = 0.5_wp * ( wfreq(ji,jj)*wfreq(ji,jj) + wfreq(ji+1,jj)*wfreq(ji+1,jj) ) / grav 137 zk_v(ji,jj) = 0.5_wp * ( wfreq(ji,jj)*wfreq(ji,jj) + wfreq(ji,jj+1)*wfreq(ji,jj+1) ) / grav 138 ! 139 zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) ) 140 zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) ) 141 END DO 142 END DO 143 ENDIF 127 144 ! 128 145 ! !== horizontal Stokes Drift 3D velocity ==! 129 DO jj = 1, jpjm1 ! exp. wave number & Stokes drift velocity at u- & v-points 130 DO ji = 1, jpim1 131 zstokes_psi_u_top(ji,jj) = 0._wp 132 zstokes_psi_v_top(ji,jj) = 0._wp 133 END DO 134 END DO 135 136 DO jk = 1, jpkm1 137 DO jj = 2, jpjm1 138 DO ji = 2, jpim1 139 zbot_u = 0.5_wp * ( gdepw_n(ji,jj,jk+1) + gdepw_n(ji+1,jj,jk+1) ) 140 zbot_v = 0.5_wp * ( gdepw_n(ji,jj,jk+1) + gdepw_n(ji,jj+1,jk+1) ) 141 zkb_u = 2.0_wp * zk_u(ji,jj) * zbot_u ! 2k * bottom depth 142 zkb_v = 2.0_wp * zk_v(ji,jj) * zbot_v ! 2k * bottom depth 143 ! 144 zke3_u = MAX(1.e-8_wp, 2.0_wp * zk_u(ji,jj) * e3u_n(ji,jj,jk)) ! 2k * thickness 145 zke3_v = MAX(1.e-8_wp, 2.0_wp * zk_v(ji,jj) * e3v_n(ji,jj,jk)) ! 2k * thickness 146 147 ! Depth attenuation .... do u component first.. 148 zdepth=zkb_u 149 zsqrt_depth = SQRT(zdepth) 150 zexp_depth = EXP(-zdepth) 151 zstokes_psi_u_bot = 1.0_wp - zexp_depth & 152 & - z_two_thirds*(zsqrtpi*zsqrt_depth*zdepth*ERFC(zsqrt_depth) & 153 & + 1.0_wp - (1.0_wp + zdepth)*zexp_depth) 154 zda_u = (zstokes_psi_u_bot - zstokes_psi_u_top(ji,jj))/zke3_u 155 zstokes_psi_u_top(ji,jj) = zstokes_psi_u_bot 156 157 ! ... and then v component 158 zdepth=zkb_v 159 zsqrt_depth = SQRT(zdepth) 160 zexp_depth = EXP(-zdepth) 161 zstokes_psi_v_bot = 1.0_wp - zexp_depth & 162 & - z_two_thirds*(zsqrtpi*zsqrt_depth*zdepth*ERFC(zsqrt_depth) & 163 & + 1.0_wp - (1.0_wp + zdepth)*zexp_depth) 164 zda_v = (zstokes_psi_v_bot - zstokes_psi_v_top(ji,jj))/zke3_v 165 zstokes_psi_v_top(ji,jj) = zstokes_psi_v_bot 166 ! 167 usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk) 168 vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk) 169 END DO 170 END DO 171 END DO 146 IF( nn_sdrift==jp_breivik ) THEN 147 DO jk = 1, jpkm1 148 DO jj = 2, jpjm1 149 DO ji = 2, jpim1 150 zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) ) 151 zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) ) 152 ! 153 zkh_u = zk_u(ji,jj) * zdep_u ! k * depth 154 zkh_v = zk_v(ji,jj) * zdep_v 155 ! ! Depth attenuation 156 zda_u = EXP( -2.0_wp*zkh_u ) / ( 1.0_wp + 8.0_wp*zkh_u ) 157 zda_v = EXP( -2.0_wp*zkh_v ) / ( 1.0_wp + 8.0_wp*zkh_v ) 158 ! 159 usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk) 160 vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk) 161 END DO 162 END DO 163 END DO 164 ELSE IF( nn_sdrift==jp_phillips .OR. nn_sdrift==jp_peakfr ) THEN 165 DO jk = 1, jpkm1 166 DO jj = 2, jpjm1 167 DO ji = 2, jpim1 168 zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) ) 169 zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) ) 170 ! 171 zkh_u = zk_u(ji,jj) * zdep_u ! k * depth 172 zkh_v = zk_v(ji,jj) * zdep_v 173 ! ! Depth attenuation 174 zda_u = EXP( -2.0_wp*zkh_u ) - SQRT(2.0_wp*rpi*zkh_u) * ERFC(SQRT(2.0_wp*zkh_u)) 175 zda_v = EXP( -2.0_wp*zkh_v ) - SQRT(2.0_wp*rpi*zkh_v) * ERFC(SQRT(2.0_wp*zkh_v)) 176 ! 177 usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk) 178 vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk) 179 END DO 180 END DO 181 END DO 182 ENDIF 183 172 184 CALL lbc_lnk( usd(:,:,:), 'U', -1. ) 173 185 CALL lbc_lnk( vsd(:,:,:), 'V', -1. ) … … 228 240 229 241 242 SUBROUTINE sbc_wstress( ) 243 !!--------------------------------------------------------------------- 244 !! *** ROUTINE sbc_wstress *** 245 !! 246 !! ** Purpose : Updates the ocean momentum modified by waves 247 !! 248 !! ** Method : - Calculate u,v components of stress depending on stress 249 !! model 250 !! - Calculate the stress module 251 !! - The wind module is not modified by waves 252 !! ** action 253 !!--------------------------------------------------------------------- 254 INTEGER :: jj, ji ! dummy loop argument 255 ! 256 IF( ln_tauoc ) THEN 257 utau(:,:) = utau(:,:)*tauoc_wave(:,:) 258 vtau(:,:) = vtau(:,:)*tauoc_wave(:,:) 259 taum(:,:) = taum(:,:)*tauoc_wave(:,:) 260 ENDIF 261 ! 262 IF( ln_tauw ) THEN 263 DO jj = 1, jpjm1 264 DO ji = 1, jpim1 265 ! Stress components at u- & v-points 266 utau(ji,jj) = 0.5_wp * ( tauw_x(ji,jj) + tauw_x(ji+1,jj) ) 267 vtau(ji,jj) = 0.5_wp * ( tauw_y(ji,jj) + tauw_y(ji,jj+1) ) 268 ! 269 ! Stress module at t points 270 taum(ji,jj) = SQRT( tauw_x(ji,jj)*tauw_x(ji,jj) + tauw_y(ji,jj)*tauw_y(ji,jj) ) 271 END DO 272 END DO 273 CALL lbc_lnk_multi( utau(:,:), 'U', -1. , vtau(:,:), 'V', -1. , taum(:,:) , 'T', -1. ) 274 ENDIF 275 ! 276 END SUBROUTINE sbc_wstress 277 278 230 279 SUBROUTINE sbc_wave( kt ) 231 280 !!--------------------------------------------------------------------- … … 250 299 ENDIF 251 300 252 IF( ln_tauoc .AND. .NOT. cpl_ wstrf) THEN !== Wave induced stress ==!301 IF( ln_tauoc .AND. .NOT. cpl_tauoc ) THEN !== Wave induced stress ==! 253 302 CALL fld_read( kt, nn_fsbc, sf_tauoc ) ! read wave norm stress from external forcing 254 303 tauoc_wave(:,:) = sf_tauoc(1)%fnow(:,:,1) 304 ENDIF 305 306 IF( ln_tauw .AND. .NOT. cpl_tauw ) THEN !== Wave induced stress ==! 307 CALL fld_read( kt, nn_fsbc, sf_tauw ) ! read ocean stress components from external forcing (T grid) 308 tauw_x(:,:) = sf_tauw(1)%fnow(:,:,1) 309 tauw_y(:,:) = sf_tauw(2)%fnow(:,:,1) 255 310 ENDIF 256 311 … … 261 316 IF( jp_hsw > 0 ) hsw (:,:) = sf_sd(jp_hsw)%fnow(:,:,1) ! significant wave height 262 317 IF( jp_wmp > 0 ) wmp (:,:) = sf_sd(jp_wmp)%fnow(:,:,1) ! wave mean period 318 IF( jp_wfr > 0 ) wfreq(:,:) = sf_sd(jp_wfr)%fnow(:,:,1) ! Peak wave frequency 263 319 IF( jp_usd > 0 ) ut0sd(:,:) = sf_sd(jp_usd)%fnow(:,:,1) ! 2D zonal Stokes Drift at T point 264 320 IF( jp_vsd > 0 ) vt0sd(:,:) = sf_sd(jp_vsd)%fnow(:,:,1) ! 2D meridional Stokes Drift at T point … … 273 329 ! !== Computation of the 3d Stokes Drift ==! 274 330 ! 275 IF( jpfld == 4 ) CALL sbc_stokes() ! Calculate only if required fields are read 276 ! ! In coupled wave model-NEMO case the call is done after coupling 331 IF( ((nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) .AND. & 332 jp_hsw>0 .AND. jp_wmp>0 .AND. jp_usd>0 .AND. jp_vsd>0) .OR. & 333 (nn_sdrift==jp_peakfr .AND. jp_wfr>0 .AND. jp_usd>0 .AND. jp_vsd>0) ) & 334 CALL sbc_stokes() ! Calculate only if required fields are read 335 ! ! In coupled wave model-NEMO case the call is done after coupling 277 336 ! 278 337 ENDIF … … 299 358 !! 300 359 CHARACTER(len=100) :: cn_dir ! Root directory for location of drag coefficient files 301 TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) :: slf_i ! array of namelist informations on the fields to read360 TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) :: slf_i, slf_j ! array of namelist informations on the fields to read 302 361 TYPE(FLD_N) :: sn_cdg, sn_usd, sn_vsd, & 303 & sn_hsw, sn_wmp, sn_wnum, sn_tauoc ! informations about the fields to be read 304 ! 305 NAMELIST/namsbc_wave/ sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wnum, sn_tauoc 362 & sn_hsw, sn_wmp, sn_wfr, sn_wnum, & 363 & sn_tauoc, sn_tauwx, sn_tauwy ! informations about the fields to be read 364 ! 365 NAMELIST/namsbc_wave/ sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wfr, & 366 sn_wnum, sn_tauoc, sn_tauwx, sn_tauwy 306 367 !!--------------------------------------------------------------------- 307 368 ! … … 328 389 329 390 IF( ln_tauoc ) THEN 330 IF( .NOT. cpl_ wstrf) THEN391 IF( .NOT. cpl_tauoc ) THEN 331 392 ALLOCATE( sf_tauoc(1), STAT=ierror ) !* allocate and fill sf_wave with sn_tauoc 332 393 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' ) … … 339 400 ENDIF 340 401 402 IF( ln_tauw ) THEN 403 IF( .NOT. cpl_tauw ) THEN 404 ALLOCATE( sf_tauw(2), STAT=ierror ) !* allocate and fill sf_wave with sn_tauwx/y 405 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_tauw structure' ) 406 ! 407 ALLOCATE( slf_j(2) ) 408 slf_j(1) = sn_tauwx 409 slf_j(2) = sn_tauwy 410 ALLOCATE( sf_tauw(1)%fnow(jpi,jpj,1) ) 411 ALLOCATE( sf_tauw(2)%fnow(jpi,jpj,1) ) 412 IF( slf_j(1)%ln_tint ) ALLOCATE( sf_tauw(1)%fdta(jpi,jpj,1,2) ) 413 IF( slf_j(2)%ln_tint ) ALLOCATE( sf_tauw(2)%fdta(jpi,jpj,1,2) ) 414 CALL fld_fill( sf_tauw, (/ slf_j /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' ) 415 ENDIF 416 ALLOCATE( tauw_x(jpi,jpj) ) 417 ALLOCATE( tauw_y(jpi,jpj) ) 418 ENDIF 419 341 420 IF( ln_sdw ) THEN ! Find out how many fields have to be read from file if not coupled 342 421 jpfld=0 343 jp_usd=0 ; jp_vsd=0 ; jp_hsw=0 ; jp_wmp=0 422 jp_usd=0 ; jp_vsd=0 ; jp_hsw=0 ; jp_wmp=0 ; jp_wfr=0 344 423 IF( .NOT. cpl_sdrftx ) THEN 345 424 jpfld = jpfld + 1 … … 350 429 jp_vsd = jpfld 351 430 ENDIF 352 IF( .NOT. cpl_hsig ) THEN431 IF( .NOT. cpl_hsig .AND. (nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) ) THEN 353 432 jpfld = jpfld + 1 354 433 jp_hsw = jpfld 355 434 ENDIF 356 IF( .NOT. cpl_wper ) THEN435 IF( .NOT. cpl_wper .AND. (nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) ) THEN 357 436 jpfld = jpfld + 1 358 437 jp_wmp = jpfld 438 ENDIF 439 IF( .NOT. cpl_wfreq .AND. nn_sdrift==jp_peakfr ) THEN 440 jpfld = jpfld + 1 441 jp_wfr = jpfld 359 442 ENDIF 360 443 … … 366 449 IF( jp_hsw > 0 ) slf_i(jp_hsw) = sn_hsw 367 450 IF( jp_wmp > 0 ) slf_i(jp_wmp) = sn_wmp 451 IF( jp_wfr > 0 ) slf_i(jp_wfr) = sn_wfr 452 368 453 ALLOCATE( sf_sd(jpfld), STAT=ierror ) !* allocate and fill sf_sd with stokes drift 369 454 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' ) … … 378 463 ALLOCATE( usd (jpi,jpj,jpk), vsd (jpi,jpj,jpk), wsd(jpi,jpj,jpk) ) 379 464 ALLOCATE( hsw (jpi,jpj) , wmp (jpi,jpj) ) 465 ALLOCATE( wfreq(jpi,jpj) ) 380 466 ALLOCATE( ut0sd(jpi,jpj) , vt0sd(jpi,jpj) ) 381 467 ALLOCATE( div_sd(jpi,jpj) )
Note: See TracChangeset
for help on using the changeset viewer.
