- Timestamp:
- 2017-04-13T12:07:16+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
r7878 r7905 18 18 USE oce ! ocean variables 19 19 USE sbc_oce ! Surface boundary condition: ocean fields 20 USE zdf_oce, ONLY : ln_zdfqiao21 20 USE bdy_oce ! open boundary condition variables 22 21 USE domvvl ! domain: variable volume layers … … 39 38 LOGICAL, PUBLIC :: cpl_hsig = .FALSE. 40 39 LOGICAL, PUBLIC :: cpl_phioc = .FALSE. 41 LOGICAL, PUBLIC :: cpl_sdrftx = .FALSE. 42 LOGICAL, PUBLIC :: cpl_sdrfty = .FALSE. 40 LOGICAL, PUBLIC :: cpl_sdrft = .FALSE. 43 41 LOGICAL, PUBLIC :: cpl_wper = .FALSE. 44 42 LOGICAL, PUBLIC :: cpl_wfreq = .FALSE. … … 46 44 LOGICAL, PUBLIC :: cpl_tauoc = .FALSE. 47 45 LOGICAL, PUBLIC :: cpl_wdrag = .FALSE. 46 47 INTEGER :: nn_sdrift ! type of parameterization to calculate vertical Stokes drift 48 INTEGER, PARAMETER :: jp_breivik = 0 ! Breivik 2015: v_z=v_0*[exp(2*k*z)/(1-8*k*z)] 49 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))] 50 INTEGER, PARAMETER :: jp_peakph = 2 ! Phillips using the peak wave number read from wave model instead of the inverse depth scale 48 51 49 52 INTEGER :: jpfld ! number of files to read for stokes drift … … 103 106 ! 104 107 ! select parameterization for the calculation of vertical Stokes drift 105 SELECT CASE ( nn_sdrift ) 106 ! 107 CASE ( jp_breivik ) 108 ! exp. wave number at t-point 109 IF( nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips ) THEN ! (Eq. (19) in Breivick et al. (2014) ) 108 110 zfac = 2.0_wp * rpi / 16.0_wp 109 DO jj = 1, jpj ! exp. wave number at t-point (Eq. (19) in Breivick et al. (2014) )111 DO jj = 1, jpj 110 112 DO ji = 1, jpi 111 113 ! Stokes drift velocity estimated from Hs and Tmean … … 126 128 END DO 127 129 END DO 128 ! 129 ! !== horizontal Stokes Drift 3D velocity ==! 130 ELSE IF( nn_sdrift==jp_peakph ) THEN ! peak wave number calculated from the peak frequency received by the wave model 131 DO jj = 1, jpjm1 132 DO ji = 1, jpim1 133 zk_u(ji,jj) = 0.5_wp * ( wfreq(ji,jj)*wfreq(ji,jj) + wfreq(ji+1,jj)*wfreq(ji+1,jj) ) / grav 134 zk_v(ji,jj) = 0.5_wp * ( wfreq(ji,jj)*wfreq(ji,jj) + wfreq(ji,jj+1)*wfreq(ji,jj+1) ) / grav 135 ! 136 zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) ) 137 zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) ) 138 END DO 139 END DO 140 ENDIF 141 ! 142 ! !== horizontal Stokes Drift 3D velocity ==! 143 IF( nn_sdrift==jp_breivik ) THEN 130 144 DO jk = 1, jpkm1 131 145 DO jj = 2, jpjm1 … … 145 159 END DO 146 160 END DO 147 CASE ( jp_phillips ) 148 DO jj = 1, jpjm1 ! Peak wavenumber & Stokes drift velocity at u- & v-points 149 DO ji = 1, jpim1 150 zk_u(ji,jj) = 0.5_wp * ( wfreq(ji,jj)*wfreq(ji,jj) + wfreq(ji+1,jj)*wfreq(ji+1,jj) ) / grav 151 zk_v(ji,jj) = 0.5_wp * ( wfreq(ji,jj)*wfreq(ji,jj) + wfreq(ji,jj+1)*wfreq(ji,jj+1) ) / grav 152 ! 153 zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) ) 154 zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) ) 155 END DO 156 END DO 157 ! 158 ! !== horizontal Stokes Drift 3D velocity ==! 161 ELSE IF( nn_sdrift==jp_phillips .OR. nn_sdrift==jp_peakph ) THEN 159 162 DO jk = 1, jpkm1 160 163 DO jj = 2, jpjm1 … … 165 168 zkh_u = zk_u(ji,jj) * zdep_u ! k * depth 166 169 zkh_v = zk_v(ji,jj) * zdep_v 167 ! ! Depth attenuation : beta=1 for Phillips168 zda_u = EXP( -2.0_wp*zkh_u ) - 1.0*SQRT(2.0*rpi*zkh_u) * ERFC(SQRT(2.0*zkh_u))169 zda_v = EXP( -2.0_wp*zkh_v ) - 1.0*SQRT(2.0*rpi*zkh_v) * ERFC(SQRT(2.0*zkh_v))170 ! ! Depth attenuation 171 zda_u = EXP( -2.0_wp*zkh_u ) - SQRT(2.0_wp*rpi*zkh_u) * ERFC(SQRT(2.0_wp*zkh_u)) 172 zda_v = EXP( -2.0_wp*zkh_v ) - SQRT(2.0_wp*rpi*zkh_v) * ERFC(SQRT(2.0_wp*zkh_v)) 170 173 ! 171 174 usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk) … … 174 177 END DO 175 178 END DO 176 END SELECT179 ENDIF 177 180 178 181 CALL lbc_lnk( usd(:,:,:), 'U', vsd(:,:,:), 'V', -1. ) … … 265 268 CALL fld_read( kt, nn_fsbc, sf_phioc ) ! read wave to ocean energy from external forcing 266 269 rn_crban(:,:) = 29.0 * sf_phioc(1)%fnow(:,:,1) ! ! Alfa is phioc*sqrt(rau0/zrhoa) : rau0=water density, zhroa= air density 267 WHERE( rn_crban < -1000.0 ) rn_crban = 0.0 268 WHERE( rn_crban > 1000.0 ) rn_crban = 0.0 269 ENDIF 270 271 IF( ln_sdw ) THEN !== Computation of the 3d Stokes Drift ==! 270 WHERE( rn_crban > 1.e8 ) rn_crban = 0.0 !remove first mask mistmatch points, then cap values in case of low friction velocity 271 WHERE( rn_crban < 0.0 ) rn_crban = 0.0 272 WHERE( rn_crban > 1000.0 ) rn_crban = 1000.0 273 ENDIF 274 275 IF( ln_sdw .OR. ln_rough ) THEN !== Computation of the 3d Stokes Drift ==! 272 276 ! 273 277 IF( jpfld > 0 ) THEN ! Read from file only if the field is not coupled … … 285 289 IF( jp_wfr > 0 ) THEN 286 290 wfreq(:,:) = sf_sd(jp_wfr)%fnow(:,:,1) ! Peak wave frequency 287 WHERE( wfreq < 0.0 ) wfreq = 0.0 01288 WHERE( wfreq > 100.0 ) wfreq = 0.0 01291 WHERE( wfreq < 0.0 ) wfreq = 0.0 292 WHERE( wfreq > 100.0 ) wfreq = 0.0 289 293 ENDIF 290 294 IF( jp_usd > 0 ) THEN … … 299 303 ENDIF 300 304 ENDIF 301 ! 305 ENDIF 306 ! 307 IF( ln_sdw ) THEN 302 308 ! Read also wave number if needed, so that it is available in coupling routines 303 309 IF( ln_zdfqiao .AND. .NOT.cpl_wnum ) THEN … … 308 314 ! !== Computation of the 3d Stokes Drift ==! 309 315 ! 310 IF( (nn_sdrift==jp_breivik .AND. jp_hsw>0 .AND. jp_wmp>0 .AND. jp_usd>0 .AND. jp_vsd>0) .OR. & 311 (nn_sdrift==jp_phillips .AND. jp_wfr>0 .AND. jp_usd>0 .AND. jp_vsd>0) ) & 316 IF( ((nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) .AND. & 317 jp_hsw>0 .AND. jp_wmp>0 .AND. jp_usd>0 .AND. jp_vsd>0) .OR. & 318 (nn_sdrift==jp_peakph .AND. jp_wfr>0 .AND. jp_usd>0 .AND. jp_vsd>0) ) & 312 319 CALL sbc_stokes() ! Calculate only if required fields are read 313 320 ! ! In coupled wave model-NEMO case the call is done after coupling … … 341 348 & sn_tauoc ! informations about the fields to be read 342 349 ! 343 NAMELIST/namsbc_wave/ sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wfr, sn_wnum, sn_tauoc, sn_phioc 350 NAMELIST/namsbc_wave/ sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wfr, sn_wnum, sn_tauoc, sn_phioc, & 351 ln_cdgw, ln_sdw, ln_stcor, ln_phioc, ln_tauoc, ln_zdfqiao, ln_rough, & 352 nn_drag, nn_sdrift, nn_wmix 344 353 !!--------------------------------------------------------------------- 345 354 ! … … 353 362 IF(lwm) WRITE ( numond, namsbc_wave ) 354 363 ! 364 IF(lwp) THEN !* Parameter print 365 WRITE(numout,*) 366 WRITE(numout,*) 'sbc_wave_init: wave physics' 367 WRITE(numout,*) '~~~~~~~~' 368 WRITE(numout,*) ' Namelist namsbc_wave : set wave physics parameters' 369 WRITE(numout,*) ' Stokes drift corr. to vert. velocity ln_sdw = ', ln_sdw 370 WRITE(numout,*) ' vertical parametrization nn_sdrift = ', nn_sdrift 371 WRITE(numout,*) ' Stokes coriolis term ln_stcor = ', ln_stcor 372 WRITE(numout,*) ' wave modified ocean stress ln_tauoc = ', ln_tauoc 373 WRITE(numout,*) ' wave to ocean energy ln_phioc = ', ln_phioc 374 WRITE(numout,*) ' vertical mixing parametrization nn_wmix = ', nn_wmix 375 WRITE(numout,*) ' neutral drag coefficient ln_cdgw = ', ln_cdgw 376 WRITE(numout,*) ' momentum formulation nn_drag = ', nn_drag 377 WRITE(numout,*) ' wave roughness length modification ln_rough = ', ln_rough 378 WRITE(numout,*) ' Qiao vertical mixing formulation ln_zdfqiao = ', ln_zdfqiao 379 ENDIF 380 381 IF ( ln_wave ) THEN 382 ! Activated wave physics but no wave physics components activated 383 IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor .OR. ln_phioc .OR. ln_rough .OR. ln_zdfqiao) ) THEN 384 CALL ctl_warn( 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauoc=F, ln_stcor=F, ln_phioc=F ', & 385 'ln_rough=F, ln_zdfqiao=F' ) 386 ELSEIF (ln_stcor .AND. .NOT. ln_sdw) THEN 387 CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T') 388 ENDIF 389 IF ( ln_cdgw .AND. .NOT.(nn_drag==jp_ukmo .OR. nn_drag==jp_std .OR. nn_drag==jp_const .OR. nn_drag==jp_mcore) ) & 390 CALL ctl_stop( 'The chosen nn_drag for momentum calculation must be 0, 1, 2, or 3') 391 IF ( ln_cdgw .AND. ln_blk_core .AND. nn_drag==0 ) & 392 CALL ctl_stop( 'The chosen nn_drag for momentum calculation in core forcing must be 1, 2, or 3') 393 IF ( ln_cdgw .AND. ln_flx .AND. nn_drag==3 ) & 394 CALL ctl_stop( 'The chosen nn_drag for momentum calculation in direct forcing must be 0, 1, or 2') 395 IF( ln_phioc .AND. .NOT.(nn_wmix==jp_craigbanner .OR. nn_wmix==jp_janssen) ) & 396 CALL ctl_stop( 'The chosen nn_wmix for wave vertical mixing must be 0, or 1' ) 397 IF( ln_sdw .AND. .NOT.(nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips .OR. nn_sdrift==jp_peakph) ) & 398 CALL ctl_stop( 'The chosen nn_sdrift for Stokes drift vertical velocity must be 0, 1, or 2' ) 399 IF( ln_zdfqiao .AND. .NOT.ln_sdw ) & 400 CALL ctl_stop( 'Qiao vertical mixing can not be used without Stokes drift (ln_sdw)' ) 401 ELSE 402 IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor .OR. ln_phioc .OR. ln_rough .OR. ln_zdfqiao ) & 403 & CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ', & 404 & 'with drag coefficient (ln_cdgw =T) ' , & 405 & 'or Stokes Drift (ln_sdw=T) ' , & 406 & 'or Stokes-Coriolis term (ln_stcor=T)', & 407 & 'or ocean stress modification due to waves (ln_tauoc=T) ', & 408 & 'or wave to ocean energy modification (ln_phioc=T) ', & 409 & 'or wave surface roughness (ln_rough=T) ', & 410 & 'or Qiao vertical mixing formulation (ln_zdfqiao=T) ' ) 411 ENDIF 412 ! 355 413 IF( ln_cdgw ) THEN 356 414 IF( .NOT. cpl_wdrag ) THEN … … 360 418 ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1) ) 361 419 IF( sn_cdg%ln_tint ) ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) ) 362 CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave_init', ' Wave module', 'namsbc_wave' )420 CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' ) 363 421 ENDIF 364 422 ALLOCATE( cdn_wave(jpi,jpj) ) … … 372 430 ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1) ) 373 431 IF( sn_tauoc%ln_tint ) ALLOCATE( sf_tauoc(1)%fdta(jpi,jpj,1,2) ) 374 CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave_init', ' Wave module', 'namsbc_wave' )432 CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' ) 375 433 ENDIF 376 434 ALLOCATE( tauoc_wave(jpi,jpj) ) … … 384 442 ALLOCATE( sf_phioc(1)%fnow(jpi,jpj,1) ) 385 443 IF( sn_phioc%ln_tint ) ALLOCATE( sf_phioc(1)%fdta(jpi,jpj,1,2) ) 386 CALL fld_fill( sf_phioc, (/ sn_phioc /), cn_dir, 'sbc_wave_init', ' Wave module', 'namsbc_wave' )444 CALL fld_fill( sf_phioc, (/ sn_phioc /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' ) 387 445 ENDIF 388 446 ALLOCATE( rn_crban(jpi,jpj) ) 389 447 ENDIF 390 448 391 IF( ln_sdw ) THEN ! Find out how many fields have to be read from file if not coupled 392 jpfld=0 393 jp_usd=0 ; jp_vsd=0 ; jp_hsw=0 ; jp_wmp=0 ; jp_wfr=0 394 IF( .NOT. cpl_sdrftx ) THEN 449 ! Find out how many fields have to be read from file if not coupled 450 jpfld=0 451 jp_usd=0 ; jp_vsd=0 ; jp_hsw=0 ; jp_wmp=0 ; jp_wfr=0 452 IF( ln_sdw ) THEN 453 IF( .NOT. cpl_sdrft ) THEN 395 454 jpfld = jpfld + 1 396 455 jp_usd = jpfld 397 ENDIF398 IF( .NOT. cpl_sdrfty ) THEN399 456 jpfld = jpfld + 1 400 457 jp_vsd = jpfld 401 458 ENDIF 402 IF( .NOT. cpl_hsig ) THEN459 IF( .NOT. cpl_hsig .AND. (nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) ) THEN 403 460 jpfld = jpfld + 1 404 461 jp_hsw = jpfld 405 462 ENDIF 406 IF( .NOT. cpl_wper ) THEN463 IF( .NOT. cpl_wper .AND. (nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) ) THEN 407 464 jpfld = jpfld + 1 408 465 jp_wmp = jpfld 409 466 ENDIF 410 IF( .NOT. cpl_wfreq ) THEN467 IF( .NOT. cpl_wfreq .AND. nn_sdrift==jp_peakph ) THEN 411 468 jpfld = jpfld + 1 412 469 jp_wfr = jpfld 413 470 ENDIF 414 415 ! Read from file only the non-coupled fields 416 IF( jpfld > 0 ) THEN 417 ALLOCATE( slf_i(jpfld) ) 418 IF( jp_usd > 0 ) slf_i(jp_usd) = sn_usd 419 IF( jp_vsd > 0 ) slf_i(jp_vsd) = sn_vsd 420 IF( jp_hsw > 0 ) slf_i(jp_hsw) = sn_hsw 421 IF( jp_wmp > 0 ) slf_i(jp_wmp) = sn_wmp 422 IF( jp_wfr > 0 ) slf_i(jp_wfr) = sn_wfr 423 ALLOCATE( sf_sd(jpfld), STAT=ierror ) !* allocate and fill sf_sd with stokes drift 424 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_sd structure' ) 425 ! 426 DO ifpr= 1, jpfld 427 ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) ) 428 IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) ) 429 END DO 430 ! 431 CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave_init', 'Wave module ', 'namsbc_wave' ) 432 ENDIF 471 ENDIF 472 473 IF( ln_rough .AND. .NOT. cpl_hsig .AND. jp_hsw==0 ) THEN 474 jpfld = jpfld + 1 475 jp_hsw = jpfld 476 ENDIF 477 478 ! Read from file only the non-coupled fields 479 IF( jpfld > 0 ) THEN 480 ALLOCATE( slf_i(jpfld) ) 481 IF( jp_usd > 0 ) slf_i(jp_usd) = sn_usd 482 IF( jp_vsd > 0 ) slf_i(jp_vsd) = sn_vsd 483 IF( jp_hsw > 0 ) slf_i(jp_hsw) = sn_hsw 484 IF( jp_wmp > 0 ) slf_i(jp_wmp) = sn_wmp 485 IF( jp_wfr > 0 ) slf_i(jp_wfr) = sn_wfr 486 ALLOCATE( sf_sd(jpfld), STAT=ierror ) !* allocate and fill sf_sd with stokes drift 487 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_sd structure' ) 488 ! 489 DO ifpr= 1, jpfld 490 ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) ) 491 IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) ) 492 END DO 493 ! 494 CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' ) 495 ENDIF 496 497 IF( ln_sdw ) THEN 433 498 ALLOCATE( usd (jpi,jpj,jpk), vsd (jpi,jpj,jpk), wsd(jpi,jpj,jpk) ) 434 ALLOCATE( hsw (jpi,jpj) , wmp (jpi,jpj))499 ALLOCATE( wmp (jpi,jpj) ) 435 500 ALLOCATE( wfreq (jpi,jpj) ) 436 501 ALLOCATE( ut0sd(jpi,jpj) , vt0sd(jpi,jpj) ) … … 446 511 ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1) ) 447 512 IF( sn_wnum%ln_tint ) ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) ) 448 CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', ' Wave module', 'namsbc_wave' )513 CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'read wave input', 'namsbc_wave' ) 449 514 ENDIF 450 515 ALLOCATE( wnum(jpi,jpj) ) 516 ENDIF 517 518 IF( ln_sdw .OR. ln_rough ) THEN 519 ALLOCATE( hsw (jpi,jpj) ) 451 520 ENDIF 452 521 !
Note: See TracChangeset
for help on using the changeset viewer.
