Changeset 4644 for branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO
- Timestamp:
- 2014-05-15T15:56:53+02:00 (10 years ago)
- Location:
- branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 2 added
- 7 edited
-
DYN/dynstcor.F90 (added)
-
SBC/sbc_oce.F90 (modified) (4 diffs)
-
SBC/sbcblk_core.F90 (modified) (19 diffs)
-
SBC/sbcmod.F90 (modified) (5 diffs)
-
SBC/sbcwave.F90 (modified) (3 diffs)
-
SBC/sbcwave_ecmwf.F90 (added)
-
ZDF/zdftke.F90 (modified) (9 diffs)
-
step.F90 (modified) (1 diff)
-
step_oce.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
r4306 r4644 53 53 LOGICAL , PUBLIC :: ln_sdw !: true if 3d stokes drift from wave model 54 54 ! 55 LOGICAL , PUBLIC :: ln_wavetke = .FALSE. !: true if wave parameters are read 56 LOGICAL , PUBLIC :: ln_stcor = .FALSE. !: true if Stokes-Coriolis forcing is included 57 LOGICAL , PUBLIC :: ln_tauoc = .FALSE. !: true if wave-modified water-side stress is used 58 55 59 LOGICAL , PUBLIC :: ln_icebergs !: Icebergs 56 60 ! … … 69 73 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vtau , vtau_b !: sea surface j-stress (ocean referential) [N/m2] 70 74 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: taum !: module of sea surface stress (at T-point) [N/m2] 75 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: cdn_wave !: wave model drag coefficient [N/m2] 76 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wspd_wavepar !: wave model 10-m neutral wind [m/s] CCC 77 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tauoc_wavepar !: normalized stress into the ocean from wave model CCC 71 78 !! wndm is used onmpute surface gases exchanges in ice-free ocean or leads 72 79 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wndm !: wind speed module at T-point (=|U10m-Uoce|) [m/s] … … 131 138 ALLOCATE( rnf (jpi,jpj) , sbc_tsc (jpi,jpj,jpts) , qsr_hc (jpi,jpj,jpk) , & 132 139 & rnf_b(jpi,jpj) , sbc_tsc_b(jpi,jpj,jpts) , qsr_hc_b(jpi,jpj,jpk) , STAT=ierr(3) ) 140 ! Initialize these since it may not done elsewhere in the code. 141 rnf (:,:) = 0.0_wp 142 sbc_tsc (:,:,:) = 0.0_wp 143 qsr_hc (:,:,:) = 0.0_wp 144 rnf_b (:,:) = 0.0_wp 145 sbc_tsc_b(:,:,:) = 0.0_wp 146 qsr_hc_b (:,:,:) = 0.0_wp 133 147 ! 134 148 ALLOCATE( tprecip(jpi,jpj) , sprecip(jpi,jpj) , fr_i(jpi,jpj) , & … … 138 152 & ssu_m (jpi,jpj) , sst_m(jpi,jpj) , & 139 153 & ssv_m (jpi,jpj) , sss_m (jpi,jpj), ssh_m(jpi,jpj) , STAT=ierr(4) ) 154 ALLOCATE( cdn_wave(jpi,jpj) ) 155 ALLOCATE( wspd_wavepar(jpi,jpj) ) 140 156 ! 141 157 #if defined key_vvl -
branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
r4624 r4644 38 38 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 39 39 USE prtctl ! Print control 40 USE sbcwave,ONLY : cdn_wave !wave module40 !USE sbcwave,ONLY : cdn_wave !wave module ! moved to sbc_oce 41 41 #if defined key_lim3 || defined key_cice 42 42 USE sbc_ice ! Surface boundary condition: ice fields … … 64 64 65 65 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf ! structure of input fields (file informations, fields read) 66 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: Cd_n10 ! 10m neutral drag coefficient for output 66 67 67 68 ! !!! CORE bulk parameters … … 83 84 REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements 84 85 REAL(wp) :: rn_zu ! z(u) : height of wind measurements 86 LOGICAL :: ln_cdec = .FALSE. ! logical flag for using ec neutral wind drag coef 87 ! References : P. Janssen, 2008, ECMWF Workshop on Ocean-Atmosphere Interactions, 10-12 November p47-60 85 88 86 89 !! * Substitutions … … 137 140 TYPE(FLD_N) :: sn_tdif ! " " 138 141 NAMELIST/namsbc_core/ cn_dir , ln_2m , ln_taudif, rn_pfac, rn_efac, rn_vfac, & 139 & sn_wndi, sn_wndj, sn_humi , sn_qsr , &140 & sn_qlw , sn_tair, sn_prec , sn_snow, &142 & sn_wndi, sn_wndj, sn_humi , sn_qsr , ln_cdec, & 143 & sn_qlw , sn_tair, sn_prec , sn_snow, & 141 144 & sn_tdif, rn_zqt , ln_bulk2z, rn_zu 142 145 !!--------------------------------------------------------------------- … … 163 166 sn_qsr%ln_tint = .false. 164 167 ENDIF 168 IF( ln_cdec .AND. lwp ) WRITE(numout,*)'Using Hans Hersbach formula for drag.' 165 169 ! ! store namelist information in an array 166 170 slf_i(jp_wndi) = sn_wndi ; slf_i(jp_wndj) = sn_wndj … … 182 186 CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_core', 'flux formulation for ocean surface boundary condition', 'namsbc_core' ) 183 187 ! 188 ! Drag coefficent so we can write to disk 189 ALLOCATE(Cd_n10(jpi,jpj)) 190 184 191 sfx(:,:) = 0._wp ! salt flux; zero unless ice is present (computed in limsbc(_2).F90) 185 192 ! … … 239 246 INTEGER :: ji, jj ! dummy loop indices 240 247 REAL(wp) :: zcoef_qsatw, zztmp ! local variable 248 REAL(wp) :: ztheta ! local variable, wind direction 241 249 REAL(wp), DIMENSION(:,:), POINTER :: zwnd_i, zwnd_j ! wind speed components at T-point 242 250 REAL(wp), DIMENSION(:,:), POINTER :: zqsatw ! specific humidity at pst … … 283 291 !CDIR COLLAPSE 284 292 #endif 285 DO jj = 2, jpjm1 286 DO ji = fs_2, fs_jpim1 ! vect. opt. 287 zwnd_i(ji,jj) = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj ) + pu(ji,jj) ) ) 288 zwnd_j(ji,jj) = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji ,jj-1) + pv(ji,jj) ) ) 293 ! Wind vector at T points 294 IF ( ln_cdgw ) THEN 295 ! Use neutral 10-m wind speed from wave model if available 296 DO jj = 2, jpjm1 297 DO ji = fs_2, fs_jpim1 ! vect. opt. 298 ! Local wind direction from non-neutral 10-m wind vector (relative to grid and no current) 299 ! since we do not have the directional information from the wave model 300 ztheta = ATAN2(sf(jp_wndi)%fnow(ji,jj,1), sf(jp_wndj)%fnow(ji,jj,1)) 301 ! Wind vector magnitude from 10-m neutral wind speed from wave model 302 zwnd_i(ji,jj) = wspd_wavepar(ji,jj) * SIN(ztheta) 303 zwnd_j(ji,jj) = wspd_wavepar(ji,jj) * COS(ztheta) 304 ! Correct for surface current, 0.0 <= rn_vfac <= 1.0 305 zwnd_i(ji,jj) = (zwnd_i(ji,jj) - 0.5*rn_vfac*(pu(ji-1,jj ) + pu(ji,jj))) 306 zwnd_j(ji,jj) = (zwnd_j(ji,jj) - 0.5*rn_vfac*(pv(ji ,jj-1) + pv(ji,jj))) 307 END DO 289 308 END DO 290 END DO 309 ELSE 310 DO jj = 2, jpjm1 311 DO ji = fs_2, fs_jpim1 ! vect. opt. 312 ! Correct for surface current, 0.0 <= rn_vfac <= 1.0 313 zwnd_i(ji,jj) = ( sf(jp_wndi)%fnow(ji,jj,1) - 0.5*rn_vfac*( pu(ji-1,jj ) + pu(ji,jj) ) ) 314 zwnd_j(ji,jj) = ( sf(jp_wndj)%fnow(ji,jj,1) - 0.5*rn_vfac*( pv(ji ,jj-1) + pv(ji,jj) ) ) 315 END DO 316 END DO 317 ENDIF 291 318 CALL lbc_lnk( zwnd_i(:,:) , 'T', -1. ) 292 319 CALL lbc_lnk( zwnd_j(:,:) , 'T', -1. ) … … 353 380 354 381 ! ... tau module, i and j component 355 DO jj = 1, jpj 356 DO ji = 1, jpi 357 zztmp = rhoa * wndm(ji,jj) * Cd(ji,jj) 358 taum (ji,jj) = zztmp * wndm (ji,jj) 359 zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 360 zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 361 END DO 362 END DO 382 ! Use cdn_wave instead of Cd if we have wave parameters 383 IF ( ln_cdgw ) THEN 384 DO jj = 1, jpj 385 DO ji = 1, jpi 386 !zztmp = rhoa * wspd_wavepar(ji,jj) * Cd_n10(ji,jj) 387 zztmp = rhoa * wndm(ji,jj) * cdn_wave(ji,jj) 388 !taum(ji,jj) = tauoc_wavepar(ji,jj) 389 !taum(ji,jj) = zztmp * wspd_wavepar(ji,jj) 390 taum (ji,jj) = zztmp * wndm (ji,jj) 391 zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 392 zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 393 ENDDO 394 ENDDO 395 ELSE 396 DO jj = 1, jpj 397 DO ji = 1, jpi 398 zztmp = rhoa * wndm(ji,jj) * Cd(ji,jj) 399 taum (ji,jj) = zztmp * wndm (ji,jj) 400 zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 401 zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 402 ENDDO 403 ENDDO 404 ENDIF 363 405 364 406 ! ... add the HF tau contribution to the wind stress module? … … 379 421 CALL lbc_lnk( utau(:,:), 'U', -1. ) 380 422 CALL lbc_lnk( vtau(:,:), 'V', -1. ) 423 424 IF (ln_tauoc) THEN 425 utau(:,:) = utau(:,:)*tauoc_wavepar(:,:) 426 vtau(:,:) = vtau(:,:)*tauoc_wavepar(:,:) 427 taum(:,:) = taum(:,:)*tauoc_wavepar(:,:) 428 ENDIF 381 429 382 430 ! Turbulent fluxes over ocean … … 755 803 REAL(wp), DIMENSION(:,:), POINTER :: dT ! air/sea temperature differeence [K] 756 804 REAL(wp), DIMENSION(:,:), POINTER :: dq ! air/sea humidity difference [K] 757 REAL(wp), DIMENSION(:,:), POINTER :: Cd_n10 ! 10m neutral drag coefficient758 805 REAL(wp), DIMENSION(:,:), POINTER :: Ce_n10 ! 10m neutral latent coefficient 759 806 REAL(wp), DIMENSION(:,:), POINTER :: Ch_n10 ! 10m neutral sensible coefficient … … 775 822 ! 776 823 CALL wrk_alloc( jpi,jpj, stab ) ! integer 777 CALL wrk_alloc( jpi,jpj, dU10, dT, dq, C d_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )824 CALL wrk_alloc( jpi,jpj, dU10, dT, dq, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L ) 778 825 CALL wrk_alloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta, U_n10, xlogt, xct, zpsi_h, zpsi_m ) 779 826 780 827 !! * Start 781 828 !! Air/sea differences 782 dU10 = max(0.5, dU) ! we don't want to fall under 0.5 m/s829 dU10 = MAX(0.5, dU) ! we don't want to fall under 0.5 m/s 783 830 dT = T_a - sst ! assuming that T_a is allready the potential temp. at zzu 784 831 dq = q_a - q_sat … … 788 835 !! 789 836 !! Neutral Drag Coefficient 790 stab = 0.5 + sign(0.5,dT) ! stable : stab = 1 ; unstable : stab = 0 791 IF ( ln_cdgw ) THEN 792 cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1) 793 Cd_n10(:,:) = cdn_wave 837 stab = 0.5 + SIGN(0.5,dT) ! stable : stab = 1 ; unstable : stab = 0 838 839 ! Use drag coefficient from a wave model ... 840 IF (ln_cdgw) THEN 841 cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1) 842 Cd_n10(:,:) = cdn_wave(:,:) 794 843 ELSE 795 Cd_n10 = 1.e-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 ) ! L & Y eq. (6a) 844 ! ... or choose a drag law 845 IF (ln_cdec) THEN 846 Cd_n10 = 1E-3 * ( 1.03 + 0.04 * dU10**1.48 ) / (dU10**0.21) ! PJ (6) 847 ELSE 848 Cd_n10 = 1E-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 ) ! L & Y eq. (6a) 849 ENDIF 796 850 ENDIF 797 851 sqrt_Cd_n10 = sqrt(Cd_n10) … … 824 878 825 879 !! Updating the neutral 10m transfer coefficients : 826 Cd_n10 = 1.e-3 * (2.7/U_n10 + 0.142 + U_n10/13.09) ! L & Y eq. (6a) 880 IF (ln_cdec) THEN 881 Cd_n10 = 1E-3 * ( 1.03 + 0.04 * U_n10**1.48 ) / (U_n10**0.21) ! PJ (6) 882 ELSE 883 Cd_n10 = 1E-3 * ( 2.7/U_n10 + 0.142 + U_n10/13.09 ) ! L & Y eq. (6a) 884 ENDIF 827 885 sqrt_Cd_n10 = sqrt(Cd_n10) 828 886 Ce_n10 = 1.e-3 * (34.6 * sqrt_Cd_n10) ! L & Y eq. (6b) … … 847 905 !! 848 906 CALL wrk_dealloc( jpi,jpj, stab ) ! integer 849 CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, C d_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )907 CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L ) 850 908 CALL wrk_dealloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta, U_n10, xlogt, xct, zpsi_h, zpsi_m ) 851 909 ! … … 891 949 REAL(wp), DIMENSION(:,:), POINTER :: dT ! air/sea temperature differeence [K] 892 950 REAL(wp), DIMENSION(:,:), POINTER :: dq ! air/sea humidity difference [K] 893 REAL(wp), DIMENSION(:,:), POINTER :: Cd_n10 ! 10m neutral drag coefficient894 951 REAL(wp), DIMENSION(:,:), POINTER :: Ce_n10 ! 10m neutral latent coefficient 895 952 REAL(wp), DIMENSION(:,:), POINTER :: Ch_n10 ! 10m neutral sensible coefficient … … 911 968 IF( nn_timing == 1 ) CALL timing_start('TURB_CORE_2Z') 912 969 ! 913 CALL wrk_alloc( jpi,jpj, dU10, dT, dq, C d_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )970 CALL wrk_alloc( jpi,jpj, dU10, dT, dq, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L ) 914 971 CALL wrk_alloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta_u, zeta_t, U_n10 ) 915 972 CALL wrk_alloc( jpi,jpj, xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m ) … … 923 980 !! Neutral Drag Coefficient : 924 981 stab = 0.5 + sign(0.5,dT) ! stab = 1 if dT > 0 -> STABLE 925 IF( ln_cdgw ) THEN 926 cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1) 927 Cd_n10(:,:) = cdn_wave 982 983 ! Use drag coefficient from a wave model ... 984 IF (ln_cdgw) THEN 985 cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1) 986 Cd_n10(:,:) = cdn_wave 928 987 ELSE 929 Cd_n10 = 1.e-3*( 2.7/dU10 + 0.142 + dU10/13.09 ) 988 ! ... or choose a drag law 989 IF (ln_cdec) THEN 990 Cd_n10 = 1E-3 * ( 1.03 + 0.04 * dU10**1.48 ) / (dU10**0.21) ! PJ (6) 991 ELSE 992 Cd_n10 = 1E-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 ) ! L & Y eq. (6a) 993 ENDIF 930 994 ENDIF 931 995 sqrt_Cd_n10 = sqrt(Cd_n10) … … 996 1060 END DO 997 1061 !! 998 CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, C d_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )1062 CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L ) 999 1063 CALL wrk_dealloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta_u, zeta_t, U_n10 ) 1000 1064 CALL wrk_dealloc( jpi,jpj, xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m ) -
branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
r4624 r4644 50 50 USE timing ! Timing 51 51 USE sbcwave ! Wave module 52 USE sbcwave_ecmwf ! ECMWF wave module 52 53 53 54 IMPLICIT NONE … … 82 83 INTEGER :: icpt ! local integer 83 84 !! 84 NAMELIST/namsbc/ nn_fsbc , ln_ana , ln_flx, ln_blk_clio, ln_blk_core, ln_cpl, & 85 & ln_blk_mfs, ln_apr_dyn, nn_ice, nn_ice_embd, ln_dm2dc , ln_rnf, & 86 & ln_ssr , nn_fwb , ln_cdgw , ln_wave , ln_sdw, nn_lsm, cn_iceflx 85 NAMELIST/namsbc/ nn_fsbc , ln_ana , ln_flx , ln_blk_clio, ln_blk_core, ln_cpl, & 86 & ln_blk_mfs, ln_apr_dyn, nn_ice , nn_ice_embd, ln_dm2dc , ln_rnf, & 87 & ln_ssr , nn_fwb , ln_cdgw , ln_wave , ln_sdw , nn_lsm, & 88 & cn_iceflx , ln_stcor , ln_tauoc, ln_wavetke 87 89 INTEGER :: ios 88 90 !!---------------------------------------------------------------------- … … 135 137 WRITE(numout,*) ' closed sea (=0/1) (set in namdom) nn_closea = ', nn_closea 136 138 WRITE(numout,*) ' n. of iterations if land-sea-mask applied nn_lsm = ', nn_lsm 139 WRITE(numout,*) ' Drag coefficient from wave model ln_cdgw = ', ln_cdgw 140 WRITE(numout,*) ' Stokes-Coriolis forcing from wave model ln_stcor = ', ln_stcor 141 WRITE(numout,*) ' TKE wave breaking BC from wave model ln_wavetke = ', ln_wavetke 142 WRITE(numout,*) ' Wave-modified stress from wave model ln_tauoc = ', ln_tauoc 137 143 ENDIF 138 144 … … 230 236 ELSE 231 237 IF ( ln_cdgw .OR. ln_sdw ) & 232 & CALL ctl_ stop('Not Activated Wave Module (ln_wave=F) but &238 & CALL ctl_warn('Not Activated Wave Module (ln_wave=F) but & 233 239 & asked coupling with drag coefficient (ln_cdgw =T) or Stokes drift (ln_sdw=T) ') 234 240 ENDIF … … 314 320 CALL sbc_ssm( kt ) ! ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m) 315 321 ! ! averaged over nf_sbc time-step 322 323 IF (ln_cdgw) CALL sbc_wave_ecmwf( kt ) 324 ! Read wave parameters if Stokes-Coriolis forcing OR wave parameters for TKE are needed 325 ! Also read it if drag coefficient to ensure that we use the neutral 10m from WAM. 326 IF (ln_stcor .OR. ln_wavetke .OR. ln_tauoc .OR. ln_cdgw) THEN 327 CALL sbc_wavepar( kt ) 328 ENDIF 316 329 317 330 IF (ln_wave) CALL sbc_wave( kt ) -
branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
r4624 r4644 31 31 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_cd ! structure of input fields (file informations, fields read) Drag Coefficient 32 32 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_sd ! structure of input fields (file informations, fields read) Stokes Drift 33 REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:) :: cdn_wave34 33 REAL(wp),ALLOCATABLE,DIMENSION (:,:) :: usd2d,vsd2d,uwavenum,vwavenum 35 34 REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:,:) :: usd3d,vsd3d,wsd3d … … 61 60 USE divcur 62 61 USE wrk_nemo 62 USE sbc_oce, ONLY : cdn_wave 63 63 #if defined key_bdy 64 64 USE bdy_oce, ONLY : bdytmask … … 100 100 IF( sn_cdg%ln_tint ) ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) ) 101 101 CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 102 ALLOCATE( cdn_wave(jpi,jpj) )103 cdn_wave(:,:) = 0.0104 102 ENDIF 105 103 IF ( ln_sdw ) THEN -
branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
r4624 r4644 26 26 !! ! + cleaning of the parameters + bugs correction 27 27 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 28 !! ! 2012-12 (oyvind.breivik@ecwmf.int) Changed to wave breaking source term 28 29 !!---------------------------------------------------------------------- 29 30 #if defined key_zdftke || defined key_esopa … … 42 43 USE domvvl ! domain: variable volume layer 43 44 USE sbc_oce ! surface boundary condition: ocean 45 USE sbcwave ! wave parameters 46 USE sbcwave_ecmwf ! ECMWF wave parameters 44 47 USE zdf_oce ! vertical physics: ocean variables 45 48 USE zdfmxl ! vertical physics: mixed layer … … 78 81 LOGICAL :: ln_lc ! Langmuir cells (LC) as a source term of TKE or not 79 82 REAL(wp) :: rn_lc ! coef to compute vertical velocity of Langmuir cells 83 REAL(wp) :: rn_alpavg = 100.0_wp ! Average alpha (normalized, water-side energy flux, only used in 84 ! case we use the source term formulation with no wave model information 85 ! NB! The option of running with source term and no wave model has 86 ! not been properly tested yet but is included for consistency (Oyvind Breivik, 2013-07-16) 87 REAL(wp) :: rn_deptmaxtkew = 50.0_wp ! maximum depth [m] to be affected by TKE from breaking waves 80 88 81 89 REAL(wp) :: ri_cri ! critic Richardson number (deduced from rn_ediff and rn_ediss values) … … 85 93 86 94 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: en !: now turbulent kinetic energy [m2/s2] 95 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e0 !: top level turbulent kinetic energy [m2/s2] from waves 96 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dfac !: wave depth scaling factor 87 97 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: htau ! depth of tke penetration (nn_htau) 88 98 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dissl ! now mixing lenght of dissipation … … 109 119 !! *** FUNCTION zdf_tke_alloc *** 110 120 !!---------------------------------------------------------------------- 121 INTEGER :: iexalloc 111 122 ALLOCATE( & 112 123 #if defined key_c1d … … 118 129 & avmu_k(jpi,jpj,jpk) , avmv_k(jpi,jpj,jpk), STAT= zdf_tke_alloc ) 119 130 ! 131 IF ( ln_wavetke ) THEN 132 ALLOCATE( e0(jpi,jpj) , dfac(jpi,jpj), STAT= iexalloc ) 133 ! Setting arrays to zero initially. 134 e0(:,:) = 0.0_wp 135 dfac(:,:) = 0.0_wp 136 ! 137 zdf_tke_alloc = zdf_tke_alloc + iexalloc 138 ENDIF 120 139 IF( lk_mpp ) CALL mpp_sum ( zdf_tke_alloc ) 121 140 IF( zdf_tke_alloc /= 0 ) CALL ctl_warn('zdf_tke_alloc: failed to allocate arrays') … … 219 238 REAL(wp) :: zus , zwlc , zind ! - - 220 239 REAL(wp) :: zzd_up, zzd_lw ! - - 240 REAL(wp) :: lambda ! - - 241 REAL(wp) :: zphio_flux ! - - 242 REAL(wp) :: sbr, z0, fb ! - - 243 221 244 !!bfr REAL(wp) :: zebot ! - - 222 245 INTEGER , POINTER, DIMENSION(:,: ) :: imlc … … 239 262 ! ! Surface boundary condition on tke 240 263 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 241 DO jj = 2, jpjm1 ! en(1) = rn_ebb taum / rau0 (min value rn_emin0) 242 DO ji = fs_2, fs_jpim1 ! vector opt. 243 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 244 END DO 245 END DO 246 264 ! Wave parameters read? 265 IF ( ln_wavetke ) THEN 266 ! Breaking-wave TKE injection as flux 267 DO jj = 2, jpjm1 ! en(1) = rn_ebb taum / rau0 (min value rn_emin0) 268 DO ji = fs_2, fs_jpim1 ! vector opt. 269 ! Roughness length proportional to wave height, capped at Hs = 0.1 m 270 z0 = 0.5_wp*MAX(swh_wavepar(ji,jj), 0.1_wp) 271 lambda = 2.38_wp/z0 272 ! The energy flux from the wave model in physical units 273 zphio_flux = MAX(phioc_wavepar(ji,jj), 0.0_wp) 274 ! TKE at surface (Mellor and Blumberg, 2004) 275 e0(ji,jj) = 0.5_wp*( 15.8_wp*zphio_flux/rau0 )**0.67_wp 276 ! Reducing the surface TKE by scaling with T-depth of first vertical level 277 dfac(ji,jj) = MAX( (2.0_wp*lambda*fsdept(ji,jj,1) )/3.0_wp, 0.00001_wp ) 278 ! TKE surface boundary condition weighted by the thickness of the first level 279 en(ji,jj,1) = MAX( rn_emin0, e0(ji,jj) * ( 1.0_wp - EXP(-dfac(ji,jj)) ) / dfac(ji,jj) ) * tmask(ji,jj,1) 280 ENDDO 281 ENDDO 282 ! If wave parameters are not read, use rn_ebb, corresponding to alpha = 100.0 283 ELSE 284 DO jj = 2, jpjm1 ! en(1) = rn_ebb taum / rau0 (min value rn_emin0) 285 DO ji = fs_2, fs_jpim1 ! vector opt. 286 ! Estimate average energy flux from std value found in rn_ebb (corresponds to alpha=100.0) 287 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 288 ENDDO 289 ENDDO 290 ENDIF 291 247 292 !!bfr - start commented area 248 293 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 443 488 ENDIF 444 489 CALL lbc_lnk( en, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) 490 IF ( ln_wavetke ) THEN 491 CALL lbc_lnk( e0, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) 492 CALL lbc_lnk( dfac, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) 493 ENDIF 445 494 ! 446 495 CALL wrk_dealloc( jpi,jpj, imlc ) ! integer -
branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/step.F90
r4624 r4644 286 286 CALL dyn_ldf( kstp ) ! lateral mixing 287 287 IF( ln_neptsimp ) CALL dyn_nept_cor( kstp ) ! add Neptune velocities (simplified) 288 IF( ln_stcor ) CALL dyn_stcor ( kstp ) ! Stokes-Coriolis forcing (ln_stcor set in SBC/sbc_oce.F90) 288 289 #if defined key_agrif 289 290 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn ! momemtum sponge -
branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
r4328 r4644 23 23 24 24 USE sbcmod ! surface boundary condition (sbc routine) 25 USE sbc_oce , ONLY: ln_stcor ! seems to be necessary but should be included in sbcmod???? 25 26 USE sbcrnf ! surface boundary condition: runoff variables 26 27 USE sbccpl ! surface boundary condition: coupled formulation (call send at end of step) … … 51 52 USE dynspg ! surface pressure gradient (dyn_spg routine) 52 53 USE dynnept ! simp. form of Neptune effect(dyn_nept_cor routine) 54 USE dynstcor ! Stokes-Coriolis forcing (dyn_stcor routine) 53 55 54 56 USE dynnxt ! time-stepping (dyn_nxt routine)
Note: See TracChangeset
for help on using the changeset viewer.
