Changeset 7809 for branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO
- Timestamp:
- 2017-03-17T15:21:06+01:00 (7 years ago)
- Location:
- branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 5 edited
-
SBC/sbc_oce.F90 (modified) (1 diff)
-
SBC/sbccpl.F90 (modified) (1 diff)
-
SBC/sbcmod.F90 (modified) (3 diffs)
-
SBC/sbcwave.F90 (modified) (8 diffs)
-
ZDF/zdfgls.F90 (modified) (14 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
r7792 r7809 68 68 LOGICAL , PUBLIC :: ln_sdw !: true if 3d stokes drift from wave model 69 69 LOGICAL , PUBLIC :: ln_tauoc !: true if normalized stress from wave is used 70 LOGICAL , PUBLIC :: ln_phioc !: true if wave energy to ocean is used 70 71 LOGICAL , PUBLIC :: ln_stcor !: true if Stokes-Coriolis term is used 71 72 ! -
branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
r7797 r7809 1195 1195 ! ! ========================= ! 1196 1196 IF( srcv(jpr_tauoc)%laction .AND. ln_tauoc ) tauoc_wave(:,:) = frcv(jpr_tauoc)%z3(:,:,1) 1197 1198 ! ! ========================= ! 1199 ! ! Wave to ocean energy ! 1200 ! ! ========================= ! 1201 IF( srcv(jpr_phioc)%laction .AND. ln_phioc ) THEN 1202 rn_crban(:,:) = 29.0 * frcv(jpr_phioc)%z3(:,:,1) 1203 ENDIF 1197 1204 1198 1205 ! Fields received by SAS when OASIS coupling -
branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
r7792 r7809 90 90 & ln_ssr , nn_isf , nn_fwb, ln_cdgw , ln_wave , ln_sdw , & 91 91 & ln_tauoc , ln_stcor , nn_lsm, nn_limflx , nn_components, ln_cpl , & 92 & ln_ wavcpl , nn_drag92 & ln_phioc , ln_wavcpl , nn_drag 93 93 INTEGER :: ios 94 94 INTEGER :: ierr, ierr0, ierr1, ierr2, ierr3, jpm … … 136 136 WRITE(numout,*) ' Stokes drift corr. to vert. velocity ln_sdw = ', ln_sdw 137 137 WRITE(numout,*) ' wave modified ocean stress ln_tauoc = ', ln_tauoc 138 WRITE(numout,*) ' wave to ocean energy - wave breaking ln_phioc = ', ln_phioc 138 139 WRITE(numout,*) ' Stokes coriolis term ln_stcor = ', ln_stcor 139 140 WRITE(numout,*) ' neutral drag coefficient ln_cdgw = ', ln_cdgw … … 224 225 IF ( ln_wave ) THEN 225 226 !Activated wave module but neither drag nor stokes drift activated 226 IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor ) ) THEN227 CALL ctl_warn( 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauoc=F, ln_stcor=F ')227 IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor .OR. ln_phioc) ) THEN 228 CALL ctl_warn( 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauoc=F, ln_stcor=F, ln_phioc=F') 228 229 ELSEIF (ln_stcor .AND. .NOT. ln_sdw) THEN 229 230 CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T') 230 231 ENDIF 231 232 ELSE 232 IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor ) &233 IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor .OR. ln_phioc ) & 233 234 & CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ', & 234 235 & 'with drag coefficient (ln_cdgw =T) ' , & 235 236 & 'or Stokes Drift (ln_sdw=T) ' , & 236 237 & 'or ocean stress modification due to waves (ln_tauoc=T) ', & 237 & 'or Stokes-Coriolis term (ln_stcori=T)' ) 238 & 'or wave to ocean energy modification (ln_phioc=T) ', & 239 & 'or Stokes-Coriolis term (ln_stcor=T)' ) 238 240 ENDIF 239 241 ! ! Choice of the Surface Boudary Condition (set nsbc) -
branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
r7797 r7809 56 56 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_wn ! structure of input fields (file informations, fields read) wave number for Qiao 57 57 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_tauoc ! structure of input fields (file informations, fields read) normalized wave stress into the ocean 58 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_phioc ! structure of input fields (file informations, fields read) wave to ocean energy 58 59 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: cdn_wave !: 59 60 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: hsw, wmp, wnum !: 61 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: rn_crban !: Craig and Banner constant for surface breaking waves mixing 60 62 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: tauoc_wave !: 61 63 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: tsd2d !: … … 222 224 ENDIF 223 225 226 IF( ln_phioc .AND. .NOT. cpl_phioc ) THEN !== Wave to ocean energy ==! 227 CALL fld_read( kt, nn_fsbc, sf_phioc ) ! read wave to ocean energy from external forcing 228 rn_crban(:,:) = 29.0 * sf_phioc(1)%fnow(:,:,1) ! ! Alfa is phioc*sqrt(rw/ra)sbc_wa 229 WHERE( rn_crban < 10.0 ) rn_crban = 10.0 230 WHERE( rn_crban > 300.0 ) rn_crban = 300.0 231 ENDIF 232 224 233 IF( ln_sdw ) THEN !== Computation of the 3d Stokes Drift ==! 225 234 ! … … 267 276 CHARACTER(len=100) :: cn_dir ! Root directory for location of drag coefficient files 268 277 TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) :: slf_i ! array of namelist informations on the fields to read 269 TYPE(FLD_N) :: sn_cdg, sn_usd, sn_vsd, &278 TYPE(FLD_N) :: sn_cdg, sn_usd, sn_vsd, sn_phioc, & 270 279 & sn_hsw, sn_wmp, sn_wnum, sn_tauoc ! informations about the fields to be read 271 280 ! 272 NAMELIST/namsbc_wave/ sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wnum, sn_tauoc 281 NAMELIST/namsbc_wave/ sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wnum, sn_tauoc, sn_phioc 273 282 !!--------------------------------------------------------------------- 274 283 ! … … 285 294 IF( .NOT. cpl_wdrag ) THEN 286 295 ALLOCATE( sf_cd(1), STAT=ierror ) !* allocate and fill sf_wave with sn_cdg 287 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_ wavestructure' )296 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_cd structure' ) 288 297 ! 289 298 ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1) ) … … 297 306 IF( .NOT. cpl_tauoc ) THEN 298 307 ALLOCATE( sf_tauoc(1), STAT=ierror ) !* allocate and fill sf_wave with sn_tauoc 299 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_ wavestructure' )308 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_tauoc structure' ) 300 309 ! 301 310 ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1) ) … … 304 313 ENDIF 305 314 ALLOCATE( tauoc_wave(jpi,jpj) ) 315 ENDIF 316 317 IF( ln_phioc ) THEN 318 IF( .NOT. cpl_phioc ) THEN 319 ALLOCATE( sf_phioc(1), STAT=ierror ) !* allocate and fill sf_wave with sn_phioc 320 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_phioc structure' ) 321 ! 322 ALLOCATE( sf_phioc(1)%fnow(jpi,jpj,1) ) 323 IF( sn_phioc%ln_tint ) ALLOCATE( sf_phioc(1)%fdta(jpi,jpj,1,2) ) 324 CALL fld_fill( sf_phioc, (/ sn_phioc /), cn_dir, 'sbc_wave_init', 'Wave module', 'namsbc_wave' ) 325 ENDIF 326 ALLOCATE( rn_crban(jpi,jpj) ) 306 327 ENDIF 307 328 … … 334 355 IF( jp_wmp > 0 ) slf_i(jp_wmp) = sn_wmp 335 356 ALLOCATE( sf_sd(jpfld), STAT=ierror ) !* allocate and fill sf_sd with stokes drift 336 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_ wavestructure' )357 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_sd structure' ) 337 358 ! 338 359 DO ifpr= 1, jpfld … … 354 375 IF( ln_zdfqiao .AND. .NOT.cpl_wnum ) THEN 355 376 ALLOCATE( sf_wn(1), STAT=ierror ) !* allocate and fill sf_wave with sn_wnum 356 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable toallocate sf_w avestructure' )377 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable toallocate sf_wn structure' ) 357 378 ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1) ) 358 379 IF( sn_wnum%ln_tint ) ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) ) -
branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
r7481 r7809 24 24 USE phycst ! physical constants 25 25 USE zdfmxl ! mixed layer 26 USE sbcwave , ONLY: hsw ! significant wave height26 USE sbcwave, ONLY: hsw,rn_crban 27 27 ! 28 28 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 48 48 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustars2 !: Squared surface velocity scale at T-points 49 49 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustarb2 !: Squared bottom velocity scale at T-points 50 51 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rsbc_tke1 52 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rsbc_tke3 53 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rsbc_psi1 50 54 51 55 ! !! ** Namelist namzdf_gls ** … … 61 65 REAL(wp) :: rn_emin ! minimum value of TKE (m2/s2) 62 66 REAL(wp) :: rn_charn ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value) 63 REAL(wp) :: rn_crban 67 REAL(wp) :: rn_crban_default ! Craig and Banner constant for surface breaking waves mixing 64 68 REAL(wp) :: rn_hsro ! Minimum surface roughness 65 REAL(wp) :: rn_frac_hs ! Fraction of wave height as surface roughness (if nn_z0_met >1)69 REAL(wp) :: rn_frac_hs ! Fraction of wave height as surface roughness (if nn_z0_met = 1) 66 70 67 71 REAL(wp) :: rcm_sf = 0.73_wp ! Shear free turbulence parameters … … 94 98 REAL(wp) :: rtrans = 0.1_wp 95 99 REAL(wp) :: rc02, rc02r, rc03, rc04 ! coefficients deduced from above parameters 96 REAL(wp) :: rsbc_tke1 , rsbc_tke2, rfact_tke! - - - -97 REAL(wp) :: rsbc_psi 1, rsbc_psi2, rfact_psi ! - - - -100 REAL(wp) :: rsbc_tke1_default, rsbc_tke2, rfact_tke ! - - - - 101 REAL(wp) :: rsbc_psi2, rsbc_psi3, rfact_psi ! - - - - 98 102 REAL(wp) :: rsbc_zs1, rsbc_zs2 ! - - - - 99 103 REAL(wp) :: rc0, rc2, rc3, rf6, rcff, rc_diff ! - - - - … … 117 121 !! *** FUNCTION zdf_gls_alloc *** 118 122 !!---------------------------------------------------------------------- 119 ALLOCATE( mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) , & 120 & ustars2(jpi,jpj) , ustarb2(jpi,jpj) , STAT= zdf_gls_alloc ) 123 ALLOCATE( mxln(jpi,jpj,jpk) , zwall(jpi,jpj,jpk) , & 124 & ustars2(jpi,jpj) , ustarb2(jpi,jpj) , & 125 & rsbc_tke1(jpi,jpj), rsbc_tke3(jpi,jpj) , & 126 & rsbc_psi1(jpi,jpj), STAT= zdf_gls_alloc ) 121 127 ! 122 128 IF( lk_mpp ) CALL mpp_sum ( zdf_gls_alloc ) … … 163 169 ustars2 = 0._wp ; ustarb2 = 0._wp ; psi = 0._wp ; zwall_psi = 0._wp 164 170 171 172 ! variable initialization 173 IF( ln_phioc ) THEN 174 rsbc_tke1(:,:) = (-rsc_tke*rn_crban(:,:)/(rcm_sf*ra_sf*rl_sf))**(2._wp/3._wp) ! k_eps = 53.Dirichlet + Wave breaking 175 rsbc_tke3(:,:) = rdt * rn_crban(:,:) ! Neumann + Wave breaking 176 rsbc_psi1(:,:) = rc0**rpp * rsbc_tke1(:,:)**rmm * rl_sf**rnn ! Dirichlet + Wave breaking 177 ENDIF 178 165 179 IF( kt /= nit000 ) THEN ! restore before value to compute tke 166 180 avt (:,:,:) = avt_k (:,:,:) … … 200 214 zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 201 215 CASE ( 3 ) ! Roughness given by the wave model (coupled or read in file) 202 zhsro(:,:) = hsw(:,:) 216 WHERE( hsw == 0._wp ) ! surface roughness length according to Charnock formula when sign. wave height 0 217 zhsro = MAX(rn_charn / grav * ustars2, rn_hsro) 218 ELSEWHERE 219 zhsro = MAX(hsw, rn_hsro) 220 END WHERE 203 221 END SELECT 204 222 … … 317 335 ! 318 336 CASE ( 0 ) ! Dirichlet case 319 ! First level 320 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 321 en(:,:,1) = MAX(en(:,:,1), rn_emin) 322 z_elem_a(:,:,1) = en(:,:,1) 323 z_elem_c(:,:,1) = 0._wp 324 z_elem_b(:,:,1) = 1._wp 325 ! 326 ! One level below 327 en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2)) & 328 & / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 329 en(:,:,2) = MAX(en(:,:,2), rn_emin ) 330 z_elem_a(:,:,2) = 0._wp 331 z_elem_c(:,:,2) = 0._wp 332 z_elem_b(:,:,2) = 1._wp 333 ! 334 ! 337 IF( ln_phioc ) THEN ! wave induced mixing case with forced/coupled fields 338 ! First level 339 en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin ) 340 z_elem_a(:,:,1) = en(:,:,1) 341 z_elem_c(:,:,1) = 0._wp 342 z_elem_b(:,:,1) = 1._wp 343 344 ! One level below 345 en(:,:,2) = MAX( rsbc_tke1(:,:) * ustars2(:,:) * ((zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin ) 346 z_elem_a(:,:,2) = 0._wp 347 z_elem_c(:,:,2) = 0._wp 348 z_elem_b(:,:,2) = 1._wp 349 ELSE ! wave induced mixing case with default values 350 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default)**(2._wp/3._wp) 351 en(:,:,1) = MAX(en(:,:,1), rn_emin) 352 z_elem_a(:,:,1) = en(:,:,1) 353 z_elem_c(:,:,1) = 0._wp 354 z_elem_b(:,:,1) = 1._wp 355 ! 356 ! One level below 357 en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default * ((zhsro(:,:)+fsdepw(:,:,2)) & 358 & / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 359 en(:,:,2) = MAX(en(:,:,2), rn_emin ) 360 z_elem_a(:,:,2) = 0._wp 361 z_elem_c(:,:,2) = 0._wp 362 z_elem_b(:,:,2) = 1._wp 363 ! 364 ! 365 ENDIF 335 366 CASE ( 1 ) ! Neumann boundary condition on d(e)/dz 336 ! 337 ! Dirichlet conditions at k=1 338 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 339 en(:,:,1) = MAX(en(:,:,1), rn_emin) 340 z_elem_a(:,:,1) = en(:,:,1) 341 z_elem_c(:,:,1) = 0._wp 342 z_elem_b(:,:,1) = 1._wp 343 ! 344 ! at k=2, set de/dz=Fw 345 !cbr 346 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 347 z_elem_a(:,:,2) = 0._wp 348 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:)) )) 349 zflxs(:,:) = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) & 350 & * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf) 351 352 en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2) 353 ! 354 ! 367 IF( ln_phioc ) THEN ! Shear free case: d(e)/dz=Fw with forced/coupled fields 368 ! Dirichlet conditions at k=1 369 en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin ) 370 z_elem_a(:,:,1) = en(:,:,1) 371 z_elem_c(:,:,1) = 0._wp 372 z_elem_b(:,:,1) = 1._wp 373 ! at k=2, set de/dz=Fw 374 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 375 z_elem_a(:,:,2) = 0._wp 376 zflxs(:,:) = rsbc_tke3(:,:) * ustars2(:,:)**1.5_wp * ((zhsro(:,:)+fsdept(:,:,1) ) / zhsro(:,:) )**(1.5*ra_sf) 377 en(:,:,2) = en(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 378 ELSE ! Shear free case: d(e)/dz=Fw with default values 379 ! Dirichlet conditions at k=1 380 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default)**(2._wp/3._wp) 381 en(:,:,1) = MAX(en(:,:,1), rn_emin) 382 z_elem_a(:,:,1) = en(:,:,1) 383 z_elem_c(:,:,1) = 0._wp 384 z_elem_b(:,:,1) = 1._wp 385 ! 386 ! at k=2, set de/dz=Fw 387 !cbr 388 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 389 z_elem_a(:,:,2) = 0._wp 390 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:)) )) 391 zflxs(:,:) = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) & 392 & * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf) 393 394 en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2) 395 ! 396 ! 397 ENDIF 355 398 END SELECT 356 399 … … 539 582 ! 540 583 CASE ( 0 ) ! Dirichlet boundary conditions 541 ! 542 ! Surface value 543 zdep(:,:) = zhsro(:,:) * rl_sf ! Cosmetic 544 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 545 z_elem_a(:,:,1) = psi(:,:,1) 546 z_elem_c(:,:,1) = 0._wp 547 z_elem_b(:,:,1) = 1._wp 548 ! 549 ! One level below 550 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdepw(:,:,2)/zhsro(:,:) ))) 551 zdep(:,:) = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:) 552 psi (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 553 z_elem_a(:,:,2) = 0._wp 554 z_elem_c(:,:,2) = 0._wp 555 z_elem_b(:,:,2) = 1._wp 556 ! 557 ! 584 IF( ln_phioc ) THEN ! Wave induced mixing case 585 ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2 586 ! balance between the production and the 587 ! dissipation terms including the wave effect 588 ! Surface value 589 zdep(:,:) = zhsro(:,:) * rl_sf ! Cosmetic 590 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 591 z_elem_a(:,:,1) = psi(:,:,1) 592 z_elem_c(:,:,1) = 0._wp 593 z_elem_b(:,:,1) = 1._wp 594 ! 595 ! One level below 596 zex1 = (rmm*ra_sf+rnn) 597 zex2 = (rmm*ra_sf) 598 zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**zex1 ) / zhsro(:,:)**zex2 599 psi (:,:,2) = rsbc_psi1(:,:) * ustars2(:,:)**rmm * zdep(:,:) * tmask(:,:,1) 600 z_elem_a(:,:,2) = 0._wp 601 z_elem_c(:,:,2) = 0._wp 602 z_elem_b(:,:,2) = 1._wp 603 ! 604 ! 605 ELSE ! Wave induced mixing case with default values 606 ! Surface value 607 zdep(:,:) = zhsro(:,:) * rl_sf ! Cosmetic 608 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 609 z_elem_a(:,:,1) = psi(:,:,1) 610 z_elem_c(:,:,1) = 0._wp 611 z_elem_b(:,:,1) = 1._wp 612 ! 613 ! One level below 614 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdepw(:,:,2)/zhsro(:,:) ))) 615 zdep(:,:) = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:) 616 psi (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 617 z_elem_a(:,:,2) = 0._wp 618 z_elem_c(:,:,2) = 0._wp 619 z_elem_b(:,:,2) = 1._wp 620 ! 621 ! 622 ENDIF 558 623 CASE ( 1 ) ! Neumann boundary condition on d(psi)/dz 559 ! 560 ! Surface value: Dirichlet 561 zdep(:,:) = zhsro(:,:) * rl_sf 562 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 563 z_elem_a(:,:,1) = psi(:,:,1) 564 z_elem_c(:,:,1) = 0._wp 565 z_elem_b(:,:,1) = 1._wp 566 ! 567 ! Neumann condition at k=2 568 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 569 z_elem_a(:,:,2) = 0._wp 570 ! 571 ! Set psi vertical flux at the surface: 572 zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 573 zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 574 zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 575 zdep(:,:) = rsbc_psi1 * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 576 & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.) 577 zflxs(:,:) = zdep(:,:) * zflxs(:,:) 578 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 579 580 ! 581 ! 624 IF( ln_phioc ) THEN ! Wave induced mixing case with forced/coupled fields 625 ! 626 zdep(:,:) = rl_sf * zhsro(:,:) 627 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 628 z_elem_a(:,:,1) = psi(:,:,1) 629 z_elem_c(:,:,1) = 0._wp 630 z_elem_b(:,:,1) = 1._wp 631 ! 632 ! Neumann condition at k=2 633 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 634 z_elem_a(:,:,2) = 0._wp 635 ! 636 ! Set psi vertical flux at the surface: 637 zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(rmm*ra_sf+rnn-1._wp) / zhsro(:,:)**(rmm*ra_sf) 638 zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) & 639 & * en(:,:,1)**rmm * zdep 640 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 641 ! 642 ELSE ! Wave induced mixing case with default values 643 ! Surface value: Dirichlet 644 zdep(:,:) = zhsro(:,:) * rl_sf 645 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 646 z_elem_a(:,:,1) = psi(:,:,1) 647 z_elem_c(:,:,1) = 0._wp 648 z_elem_b(:,:,1) = 1._wp 649 ! 650 ! Neumann condition at k=2 651 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 652 z_elem_a(:,:,2) = 0._wp 653 ! 654 ! Set psi vertical flux at the surface: 655 zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 656 zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 657 zflxs(:,:) = (rnn + rsbc_tke1_default * (rnn + rmm*ra_sf) * zdep(:,:)) * & 658 (1._wp + rsbc_tke1_default * zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 659 zdep(:,:) = rsbc_psi1(:,:) * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 660 & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.) 661 zflxs(:,:) = zdep(:,:) * zflxs(:,:) 662 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 663 ! 664 ! 665 ENDIF 582 666 END SELECT 583 667 … … 870 954 NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, & 871 955 & rn_clim_galp, ln_sigpsi, rn_hsro, & 872 & rn_crban , rn_charn, rn_frac_hs,&956 & rn_crban_default, rn_charn, rn_frac_hs,& 873 957 & nn_bc_surf, nn_bc_bot, nn_z0_met, & 874 958 & nn_stab_func, nn_clos … … 898 982 WRITE(numout,*) ' TKE Bottom boundary condition nn_bc_bot = ', nn_bc_bot 899 983 WRITE(numout,*) ' Modify psi Schmidt number (wb case) ln_sigpsi = ', ln_sigpsi 900 WRITE(numout,*) ' Craig and Banner coefficient rn_crban = ', rn_crban984 WRITE(numout,*) ' Craig and Banner coefficient (default) rn_crban = ', rn_crban_default 901 985 WRITE(numout,*) ' Charnock coefficient rn_charn = ', rn_charn 902 986 WRITE(numout,*) ' Surface roughness formula nn_z0_met = ', nn_z0_met … … 915 999 IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' ) 916 1000 IF( nn_z0_met < 0 .OR. nn_z0_met > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_z0_met is 0, 1, 2 or 3' ) 917 IF( nn_z0_met == 3 .AND. .NOT.ln_ wave ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_wave=T' )1001 IF( nn_z0_met == 3 .AND. .NOT.ln_sdw ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_sdw=T' ) 918 1002 IF( nn_stab_func < 0 .OR. nn_stab_func > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_stab_func is 0, 1, 2 and 3' ) 919 1003 IF( nn_clos < 0 .OR. nn_clos > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' ) … … 1089 1173 & - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) ) 1090 1174 1091 IF ( rn_crban==0._wp ) THEN1175 IF( .NOT. ln_phioc .AND. rn_crban_default==0._wp ) THEN 1092 1176 rl_sf = vkarmn 1093 1177 ELSE … … 1128 1212 rc03 = rc02 * rc0 1129 1213 rc04 = rc03 * rc0 1130 rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf ! Dirichlet + Wave breaking1131 rsbc_tke2 = rdt * rn_crban / rl_sf ! Neumann + Wave breaking1132 zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp )1133 rtrans = 0.2_wp / zcr ! Ad. inverse transition length between log and wave layer1214 rsbc_tke1_default = -3._wp/2._wp*rn_crban_default*ra_sf*rl_sf 1215 rsbc_tke2 = rdt * rn_crban_default / rl_sf 1216 zcr = MAX( rsmall, rsbc_tke1_default**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 1217 rtrans = 0.2_wp / zcr 1134 1218 rsbc_zs1 = rn_charn/grav ! Charnock formula for surface roughness 1135 1219 rsbc_zs2 = rn_frac_hs / 0.85_wp / grav * 665._wp ! Rascle formula for surface roughness
Note: See TracChangeset
for help on using the changeset viewer.
