Changeset 12991 for NEMO/branches
- Timestamp:
- 2020-05-29T12:58:31+02:00 (4 years ago)
- Location:
- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves
- Files:
-
- 13 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg
r12501 r12991 89 89 ! ! =2 annual global mean of e-p-r set to zero 90 90 ln_wave = .false. ! Activate coupling with wave (T => fill namsbc_wave) 91 ln_cdgw = .false. ! Neutral drag coefficient read from wave model (T => ln_wave=.true. & fill namsbc_wave)92 ln_sdw = .false. ! Read 2D Surf Stokes Drift & Computation of 3D stokes drift (T => ln_wave=.true. & fill namsbc_wave)93 nn_sdrift = 0 ! Parameterization for the calculation of 3D-Stokes drift from the surface Stokes drift94 ! ! = 0 Breivik 2015 parameterization: v_z=v_0*[exp(2*k*z)/(1-8*k*z)]95 ! ! = 1 Phillips: v_z=v_o*[exp(2*k*z)-beta*sqrt(-2*k*pi*z)*erfc(sqrt(-2*k*z))]96 ! ! = 2 Phillips as (1) but using the wave frequency from a wave model97 ln_tauwoc = .false. ! Activate ocean stress modified by external wave induced stress (T => ln_wave=.true. & fill namsbc_wave)98 ln_tauw = .false. ! Activate ocean stress components from wave model99 ln_stcor = .false. ! Activate Stokes Coriolis term (T => ln_wave=.true. & ln_sdw=.true. & fill namsbc_wave)100 91 / 101 92 !----------------------------------------------------------------------- … … 162 153 &namsbc_wave ! External fields from wave model (ln_wave=T) 163 154 !----------------------------------------------------------------------- 155 ln_sdw = .false. ! get the 2D Surf Stokes Drift & Compute the 3D stokes drift 156 ln_stcor = .false. ! add Stokes Coriolis and tracer advection terms 157 ln_cdgw = .false. ! Neutral drag coefficient read from wave model 158 ln_tauoc = .false. ! ocean stress is modified by wave induced stress 159 ln_wave_test= .false. ! Test case with constant wave fields 160 ! 161 ln_charn = .false. ! Charnock coefficient read from wave model (IFS only) 162 ln_taw = .false. ! ocean stress is modified by wave induced stress (coupled mode) 163 ln_phioc = .false. ! TKE flux from wave model 164 ln_bern_srfc= .false. ! wave induced pressure. Bernoulli head J term 165 ln_breivikFV_2016 = .false. ! breivik 2016 vertical stokes profile 166 ln_vortex_force = .false. 167 ! 168 cn_dir = './' ! root directory for the waves data location 169 !___________!_________________________!___________________!___________!_____________!________!___________!__________________!__________!_______________! 170 ! ! file name ! frequency (hours) ! variable ! time interp.! clim ! 'yearly'/ ! weights filename ! rotation ! land/sea mask ! 171 ! ! ! (if <0 months) ! name ! (logical) ! (T/F) ! 'monthly' ! ! pairing ! filename ! 172 sn_cdg = 'sdw_ecwaves_orca2' , 6. , 'drag_coeff' , .true. , .true. , 'yearly' , '' , '' , '' 173 sn_usd = 'sdw_ecwaves_orca2' , 6. , 'u_sd2d' , .true. , .true. , 'yearly' , '' , '' , '' 174 sn_vsd = 'sdw_ecwaves_orca2' , 6. , 'v_sd2d' , .true. , .true. , 'yearly' , '' , '' , '' 175 sn_hsw = 'sdw_ecwaves_orca2' , 6. , 'hs' , .true. , .true. , 'yearly' , '' , '' , '' 176 sn_wmp = 'sdw_ecwaves_orca2' , 6. , 'wmp' , .true. , .true. , 'yearly' , '' , '' , '' 177 sn_wnum = 'sdw_ecwaves_orca2' , 6. , 'wave_num' , .true. , .true. , 'yearly' , '' , '' , '' 164 178 / 165 179 !----------------------------------------------------------------------- … … 374 388 ! = 3 as = 1 applied on HF part of the stress (ln_cpl=T) 375 389 rn_eice = 0 ! below sea ice: =0 ON ; =4 OFF when ice fraction > 1/4 390 ln_mxhsw = .false. ! surface mixing length scale = F(wave height) 376 391 / 377 392 !----------------------------------------------------------------------- -
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/cfgs/SHARED/namelist_ref
r12530 r12991 230 230 ln_apr_dyn = .false. ! Patm gradient added in ocean & ice Eqs. (T => fill namsbc_apr ) 231 231 ln_wave = .false. ! Activate coupling with wave (T => fill namsbc_wave) 232 ln_cdgw = .false. ! Neutral drag coefficient read from wave model (T => ln_wave=.true. & fill namsbc_wave)233 ln_sdw = .false. ! Read 2D Surf Stokes Drift & Computation of 3D stokes drift (T => ln_wave=.true. & fill namsbc_wave)234 nn_sdrift = 0 ! Parameterization for the calculation of 3D-Stokes drift from the surface Stokes drift235 ! ! = 0 Breivik 2015 parameterization: v_z=v_0*[exp(2*k*z)/(1-8*k*z)]236 ! ! = 1 Phillips: v_z=v_o*[exp(2*k*z)-beta*sqrt(-2*k*pi*z)*erfc(sqrt(-2*k*z))]237 ! ! = 2 Phillips as (1) but using the wave frequency from a wave model238 ln_tauwoc = .false. ! Activate ocean stress modified by external wave induced stress (T => ln_wave=.true. & fill namsbc_wave)239 ln_tauw = .false. ! Activate ocean stress components from wave model240 ln_stcor = .false. ! Activate Stokes Coriolis term (T => ln_wave=.true. & ln_sdw=.true. & fill namsbc_wave)241 232 nn_lsm = 0 ! =0 land/sea mask for input fields is not applied (keep empty land/sea mask filename field) , 242 233 ! =1:n number of iterations of land/sea mask application for input fields (fill land/sea mask filename field) … … 362 353 sn_rcv_cal = 'coupled' , 'no' , '' , '' , '' 363 354 sn_rcv_co2 = 'coupled' , 'no' , '' , '' , '' 364 sn_rcv_hsig = 'none' , 'no' , '' , '' , ''365 355 sn_rcv_iceflx = 'none' , 'no' , '' , '' , '' 366 356 sn_rcv_mslp = 'none' , 'no' , '' , '' , '' 367 sn_rcv_phioc = 'none' , 'no' , '' , '' , ''368 sn_rcv_sdrfx = 'none' , 'no' , '' , '' , ''369 sn_rcv_sdrfy = 'none' , 'no' , '' , '' , ''370 sn_rcv_wper = 'none' , 'no' , '' , '' , ''371 sn_rcv_wnum = 'none' , 'no' , '' , '' , ''372 sn_rcv_wfreq = 'none' , 'no' , '' , '' , ''373 sn_rcv_wdrag = 'none' , 'no' , '' , '' , ''374 357 sn_rcv_ts_ice = 'none' , 'no' , '' , '' , '' 375 358 sn_rcv_isf = 'none' , 'no' , '' , '' , '' 376 359 sn_rcv_icb = 'none' , 'no' , '' , '' , '' 377 sn_rcv_tauwoc = 'none' , 'no' , '' , '' , '' 378 sn_rcv_tauw = 'none' , 'no' , '' , '' , '' 379 sn_rcv_wdrag = 'none' , 'no' , '' , '' , '' 360 sn_rcv_hsig = 'none' , 'no' , '' ' '' , 'T' 361 sn_rcv_phioc = 'none' , 'no' , '' , '' , 'T' 362 sn_rcv_sdrfx = 'none' , 'no' , '' , '' , 'T' 363 sn_rcv_sdrfy = 'none' , 'no' , '' ' '' , 'T' 364 sn_rcv_wper = 'none' , 'no' , '' ' '' , 'T' 365 sn_rcv_wnum = 'none' , 'no' , '' ' '' , 'T' 366 sn_rcv_wstrf = 'none' , 'no' , '' ' '' , 'T' 367 sn_rcv_wdrag = 'none' , 'no' , '' ' '' , 'T' 368 sn_rcv_charn = 'none' , 'no' , '' , '' , 'T' 369 sn_rcv_taw = 'none' , 'no' , '' , '' , 'U,V' 370 sn_rcv_bhd = 'none' , 'no' , '' ' '' , 'T' 371 sn_rcv_tusd = 'none' , 'no' , '' ' '' , 'T' 372 sn_rcv_tvsd = 'none' , 'no' , '' ' '' , 'T' 380 373 / 381 374 !----------------------------------------------------------------------- … … 555 548 &namsbc_wave ! External fields from wave model (ln_wave=T) 556 549 !----------------------------------------------------------------------- 550 ln_sdw = .false. ! get the 2D Surf Stokes Drift & Compute the 3D stokes drift 551 ln_stcor = .false. ! add Stokes Coriolis and tracer advection terms 552 ln_cdgw = .false. ! Neutral drag coefficient read from wave model 553 ln_tauoc = .false. ! ocean stress is modified by wave induced stress 554 ln_wave_test= .false. ! Test case with constant wave fields 555 ! 556 ln_charn = .false. ! Charnock coefficient read from wave model (IFS only) 557 ln_taw = .false. ! ocean stress is modified by wave induced stress (coupled mode) 558 ln_phioc = .false. ! TKE flux from wave model 559 ln_bern_srfc= .false. ! wave induced pressure. Bernoulli head J term 560 ln_breivikFV_2016 = .false. ! breivik 2016 vertical stokes profile 561 ln_vortex_force = .false. ! Vortex Force term 562 ln_stshear = .false. ! include stokes shear in EKE computation 563 ! 557 564 cn_dir = './' ! root directory for the waves data location 558 565 !___________!_________________________!___________________!___________!_____________!________!___________!__________________!__________!_______________! … … 564 571 sn_hsw = 'sdw_ecwaves_orca2' , 6. , 'hs' , .true. , .true. , 'yearly' , '' , '' , '' 565 572 sn_wmp = 'sdw_ecwaves_orca2' , 6. , 'wmp' , .true. , .true. , 'yearly' , '' , '' , '' 566 sn_wfr = 'sdw_ecwaves_orca2' , 6. , 'wfr' , .true. , .true. , 'yearly' , '' , '' , ''567 573 sn_wnum = 'sdw_ecwaves_orca2' , 6. , 'wave_num' , .true. , .true. , 'yearly' , '' , '' , '' 568 sn_tauwoc = 'sdw_ecwaves_orca2' , 6. , 'wave_stress', .true. , .true. , 'yearly' , '' , '' , ''569 sn_tauwx = 'sdw_ecwaves_orca2' , 6. , 'wave_stress', .true. , .true. , 'yearly' , '' , '' , ''570 sn_tauwy = 'sdw_ecwaves_orca2' , 6. , 'wave_stress', .true. , .true. , 'yearly' , '' , '' , ''571 574 / 572 575 !----------------------------------------------------------------------- … … 1135 1138 ln_mxl0 = .true. ! surface mixing length scale = F(wind stress) (T) or not (F) 1136 1139 rn_mxl0 = 0.04 ! surface buoyancy lenght scale minimum value 1140 ln_mxhsw = .false. ! surface mixing length scale = F(wave height) 1137 1141 ln_drg = .false. ! top/bottom friction added as boundary condition of TKE 1138 1142 ln_lc = .true. ! Langmuir cell parameterisation (Axell 2002) … … 1147 1151 ! = 1 0.5m at the equator to 30m poleward of 40 degrees 1148 1152 rn_eice = 4 ! below sea ice: =0 ON ; =4 OFF when ice fraction > 1/4 1153 nn_bc_surf = 1 ! surface condition (0/1=Dir/Neum) ! Only applicable for wave coupling (ln_cplwave=1) 1154 nn_bc_bot = 1 ! bottom condition (0/1=Dir/Neum) ! Only applicable for wave coupling (ln_cplwave=1) 1149 1155 / 1150 1156 !----------------------------------------------------------------------- -
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DYN/dynspg.F90
r12489 r12991 19 19 USE sbc_ice , ONLY : snwice_mass, snwice_mass_b 20 20 USE sbcapr ! surface boundary condition: atmospheric pressure 21 USE sbcwave, ONLY : bhd_wave 21 22 USE dynspg_exp ! surface pressure gradient (dyn_spg_exp routine) 22 23 USE dynspg_ts ! surface pressure gradient (dyn_spg_ts routine) … … 143 144 ENDIF 144 145 ! 146 IF( ln_wave .and. ln_bern_srfc ) THEN !== Add J terms: depth-independent Bernoulli head 147 DO_2D_00_00 148 spgu(ji,jj) = spgu(ji,jj) + ( bhd_wave(ji+1,jj) - bhd_wave(ji,jj) ) / e1u(ji,jj) !++ bhd_wave from wave model in m2/s2 [BHD parameters in WW3] 149 spgv(ji,jj) = spgv(ji,jj) + ( bhd_wave(ji,jj+1) - bhd_wave(ji,jj) ) / e2v(ji,jj) 150 END_2D 151 ENDIF 152 ! 145 153 DO_3D_00_00( 1, jpkm1 ) 146 154 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + spgu(ji,jj) -
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DYN/dynvor.F90
r12377 r12991 21 21 !! - ! 2018-03 (G. Madec) add two new schemes (ln_dynvor_enT and ln_dynvor_eet) 22 22 !! - ! 2018-04 (G. Madec) add pre-computed gradient for metric term calculation 23 !! 4.2 ! 2020-05 (G. Madec, E. Clementi) add vortex force trends (ln_vortex_force=T) 23 24 !!---------------------------------------------------------------------- 24 25 … … 37 38 USE trddyn ! trend manager: dynamics 38 39 USE sbcwave ! Surface Waves (add Stokes-Coriolis force) 39 USE sbc_oce , ONLY : ln_stcor! use Stoke-Coriolis force40 USE sbc_oce, ONLY : ln_stcor, ln_vortex_force ! use Stoke-Coriolis force 40 41 ! 41 42 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 119 120 ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 120 121 ! 121 ztrdu(:,:,:) = puu(:,:,:,Krhs) !* planetary vorticity trend (including Stokes-Coriolis force)122 ztrdu(:,:,:) = puu(:,:,:,Krhs) !* planetary vorticity trend 122 123 ztrdv(:,:,:) = pvv(:,:,:,Krhs) 123 124 SELECT CASE( nvor_scheme ) 124 125 CASE( np_ENS ) ; CALL vor_ens( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! enstrophy conserving scheme 125 IF( ln_stcor ) CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend126 126 CASE( np_ENE, np_MIX ) ; CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! energy conserving scheme 127 IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend128 127 CASE( np_ENT ) ; CALL vor_enT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! energy conserving scheme (T-pts) 129 IF( ln_stcor ) CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend130 128 CASE( np_EET ) ; CALL vor_eeT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! energy conserving scheme (een with e3t) 131 IF( ln_stcor ) CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend132 129 CASE( np_EEN ) ; CALL vor_een( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! energy & enstrophy scheme 133 IF( ln_stcor ) CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend134 130 END SELECT 135 131 ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:) … … 159 155 CASE( np_ENT ) !* energy conserving scheme (T-pts) 160 156 CALL vor_enT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 161 IF( ln_stcor ) CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 157 IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN 158 CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 159 ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN 160 CALL vor_enT( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force 161 ENDIF 162 162 CASE( np_EET ) !* energy conserving scheme (een scheme using e3t) 163 163 CALL vor_eeT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 164 IF( ln_stcor ) CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 164 IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN 165 CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 166 ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN 167 CALL vor_eeT( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force 168 ENDIF 165 169 CASE( np_ENE ) !* energy conserving scheme 166 170 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 167 IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 171 IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN 172 CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 173 ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN 174 CALL vor_ene( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force 175 ENDIF 168 176 CASE( np_ENS ) !* enstrophy conserving scheme 169 177 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 170 IF( ln_stcor ) CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 178 179 IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN 180 CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 181 ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN 182 CALL vor_ens( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force 183 ENDIF 171 184 CASE( np_MIX ) !* mixed ene-ens scheme 172 185 CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! relative vorticity or metric trend (ens) 173 186 CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! planetary vorticity trend (ene) 174 IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 187 IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 188 IF( ln_vortex_force ) CALL vor_ens( kt, Kmm, nrvm, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add vortex force 175 189 CASE( np_EEN ) !* energy and enstrophy conserving scheme 176 190 CALL vor_een( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 177 IF( ln_stcor ) CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 191 IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN 192 CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 193 ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN 194 CALL vor_een( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force 195 ENDIF 178 196 END SELECT 179 197 ! -
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DYN/dynzad.F90
r12377 r12991 16 16 USE trd_oce ! trends: ocean variables 17 17 USE trddyn ! trend manager: dynamics 18 USE sbcwave, ONLY: wsd ! Surface Waves (add vertical Stokes-drift) 18 19 ! 19 20 USE in_out_manager ! I/O manager … … 78 79 DO jk = 2, jpkm1 ! Vertical momentum advection at level w and u- and v- vertical 79 80 DO_2D_01_01 81 IF( ln_vortex_force ) THEN 82 zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * ( ww(ji,jj,jk) + wsd(ji,jj,jk) ) 83 ELSE 80 84 zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk) 85 ENDIF 81 86 END_2D 82 87 DO_2D_00_00 -
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/SBC/sbc_oce.F90
r12377 r12991 12 12 !! 4.0 ! 2016-06 (L. Brodeau) new unified bulk routine (based on AeroBulk) 13 13 !! 4.0 ! 2019-03 (F. Lemarié, G. Samson) add compatibility with ABL mode 14 !! 4.2 ! 2020-05 9 Madec, E. Clementi) modified wave parameters in namelist 14 15 !!---------------------------------------------------------------------- 15 16 … … 36 37 LOGICAL , PUBLIC :: ln_blk !: bulk formulation 37 38 LOGICAL , PUBLIC :: ln_abl !: Atmospheric boundary layer model 39 LOGICAL , PUBLIC :: ln_wave !: wave in the system (forced or coupled) 38 40 #if defined key_oasis3 39 41 LOGICAL , PUBLIC :: lk_oasis = .TRUE. !: OASIS used … … 56 58 ! !: = 1 global mean of e-p-r set to zero at each nn_fsbc time step 57 59 ! !: = 2 annual global mean of e-p-r set to zero 58 LOGICAL , PUBLIC :: ln_wave !: true if some coupling with wave model59 LOGICAL , PUBLIC :: ln_cdgw !: true if neutral drag coefficient from wave model60 LOGICAL , PUBLIC :: ln_sdw !: true if 3d stokes drift from wave model61 LOGICAL , PUBLIC :: ln_tauwoc !: true if normalized stress from wave is used62 LOGICAL , PUBLIC :: ln_tauw !: true if ocean stress components from wave is used63 LOGICAL , PUBLIC :: ln_stcor !: true if Stokes-Coriolis term is used64 !65 INTEGER , PUBLIC :: nn_sdrift ! type of parameterization to calculate vertical Stokes drift66 !67 60 LOGICAL , PUBLIC :: ln_icebergs !: Icebergs 68 61 ! … … 71 64 ! !!* namsbc_cpl namelist * 72 65 INTEGER , PUBLIC :: nn_cats_cpl !: Number of sea ice categories over which the coupling is carried out 73 66 ! 67 ! !!* namsbc_wave namelist * 68 LOGICAL , PUBLIC :: ln_sdw !: =T 3d stokes drift from wave model 69 LOGICAL , PUBLIC :: ln_stcor !: =T if Stokes-Coriolis and tracer advection terms are used 70 LOGICAL , PUBLIC :: ln_cdgw !: =T neutral drag coefficient from wave model 71 LOGICAL , PUBLIC :: ln_tauoc !: =T if normalized stress from wave is used 72 LOGICAL , PUBLIC :: ln_wave_test !: =T wave test case (constant Stokes drift) 73 LOGICAL , PUBLIC :: ln_charn !: =T Chranock coefficient from wave model 74 LOGICAL , PUBLIC :: ln_taw !: =T wind stress corrected by wave intake 75 LOGICAL , PUBLIC :: ln_phioc !: =T TKE surface BC from wave model 76 LOGICAL , PUBLIC :: ln_bern_srfc !: Bernoulli head, waves' inuced pressure 77 LOGICAL , PUBLIC :: ln_breivikFV_2016 !: Breivik 2016 profile 78 LOGICAL , PUBLIC :: ln_vortex_force !: vortex force activation 79 LOGICAL , PUBLIC :: ln_stshear !: Stoked Drift shear contribution in zdftke 80 ! 74 81 !!---------------------------------------------------------------------- 75 82 !! switch definition (improve readability) … … 81 88 INTEGER , PUBLIC, PARAMETER :: jp_purecpl = 5 !: Pure ocean-atmosphere Coupled formulation 82 89 INTEGER , PUBLIC, PARAMETER :: jp_none = 6 !: for OPA when doing coupling via SAS module 83 84 !!---------------------------------------------------------------------- 85 !! Stokes drift parametrization definition 86 !!---------------------------------------------------------------------- 87 INTEGER , PUBLIC, PARAMETER :: jp_breivik_2014 = 0 !: Breivik 2014: v_z=v_0*[exp(2*k*z)/(1-8*k*z)] 88 INTEGER , PUBLIC, PARAMETER :: jp_li_2017 = 1 !: Li et al 2017: Stokes drift based on Phillips spectrum (Breivik 2016) 89 ! with depth averaged profile 90 INTEGER , PUBLIC, PARAMETER :: jp_peakfr = 2 !: Li et al 2017: using the peak wave number read from wave model instead 91 ! of the inverse depth scale 92 LOGICAL , PUBLIC :: ll_st_bv2014 = .FALSE. ! logical indicator, .true. if Breivik 2014 parameterisation is active. 93 LOGICAL , PUBLIC :: ll_st_li2017 = .FALSE. ! logical indicator, .true. if Li 2017 parameterisation is active. 94 LOGICAL , PUBLIC :: ll_st_bv_li = .FALSE. ! logical indicator, .true. if either Breivik or Li parameterisation is active. 95 LOGICAL , PUBLIC :: ll_st_peakfr = .FALSE. ! logical indicator, .true. if using Li 2017 with peak wave number 96 90 ! 97 91 !!---------------------------------------------------------------------- 98 92 !! component definition -
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/SBC/sbcblk.F90
r12629 r12991 284 284 END DO 285 285 ! 286 IF( ln_wave ) THEN287 !Activated wave module but neither drag nor stokes drift activated288 IF( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) ) THEN289 CALL ctl_stop( 'STOP', 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F' )290 !drag coefficient read from wave model definable only with mfs bulk formulae and core291 ELSEIF(ln_cdgw .AND. .NOT. ln_NCAR ) THEN292 CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae')293 ELSEIF(ln_stcor .AND. .NOT. ln_sdw) THEN294 CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T')295 ENDIF296 ELSE297 IF( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) &298 & CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ', &299 & 'with drag coefficient (ln_cdgw =T) ' , &300 & 'or Stokes Drift (ln_sdw=T) ' , &301 & 'or ocean stress modification due to waves (ln_tauwoc=T) ', &302 & 'or Stokes-Coriolis term (ln_stcori=T)' )303 ENDIF304 !305 286 IF( ln_abl ) THEN ! ABL: read 3D fields for wind, temperature, humidity and pressure gradient 306 287 rn_zqt = ght_abl(2) ! set the bulk altitude to ABL first level -
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/SBC/sbcblk_algo_ecmwf.F90
r12615 r12991 17 17 !!---------------------------------------------------------------------- 18 18 !! History : 4.0 ! 2016-02 (L.Brodeau) Original code 19 !! 4.2 ! 2020-05 (G. MAdec, E. Clementi) Charnock coeff from wave model 19 20 !!---------------------------------------------------------------------- 20 21 … … 31 32 USE in_out_manager ! I/O manager 32 33 USE prtctl ! Print control 33 USE sbcwave, ONLY : cdn_wave! wave module34 USE sbcwave, ONLY : charn, ln_charn ! wave module 34 35 #if defined key_si3 || defined key_cice 35 36 USE sbc_ice ! Surface boundary condition: ice fields … … 233 234 u_star = 0.035_wp*U_blk*ztmp1/ztmp0 ! (u* = 0.035*Un10) 234 235 235 z0 = charn0*u_star*u_star/grav + 0.11_wp*znu_a/u_star 236 IF (ln_charn) THEN ! Charnock value if wave coupling 237 z0 = charn*u_star*u_star/grav + 0.11*znu_a/u_star 238 ELSE 239 z0 = charn0*u_star*u_star/grav + 0.11*znu_a/u_star 240 ENDIF 241 236 242 z0 = MIN( MAX(ABS(z0), 1.E-9) , 1._wp ) ! (prevents FPE from stupid values from masked region later on) 237 243 … … 302 308 ztmp2 = u_star*u_star 303 309 ztmp1 = znu_a/u_star 304 z0 = MIN( ABS( alpha_M*ztmp1 + charn0*ztmp2/grav ) , 0.001_wp) 310 IF (ln_charn) THEN ! Charnock value if wave coupling 311 z0 = MIN( ABS( alpha_M*ztmp1 + charn*ztmp2/grav ) , 0.001_wp) 312 ELSE 313 z0 = MIN( ABS( alpha_M*ztmp1 + charn0*ztmp2/grav ) , 0.001_wp) 314 ENDIF 305 315 z0t = MIN( ABS( alpha_H*ztmp1 ) , 0.001_wp) ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1 306 316 z0q = MIN( ABS( alpha_Q*ztmp1 ) , 0.001_wp) -
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/SBC/sbccpl.F90
r12620 r12991 8 8 !! 3.1 ! 2009_02 (G. Madec, S. Masson, E. Maisonave, A. Caubel) generic coupled interface 9 9 !! 3.4 ! 2011_11 (C. Harris) more flexibility + multi-category fields 10 !! 4.2 ! 2020-05 (G. Madec, E. Clementi) wave coupling updates 10 11 !!---------------------------------------------------------------------- 11 12 … … 102 103 INTEGER, PARAMETER :: jpr_fraqsr = 42 ! fraction of solar net radiation absorbed in the first ocean level 103 104 INTEGER, PARAMETER :: jpr_mslp = 43 ! mean sea level pressure 104 INTEGER, PARAMETER :: jpr_hsig = 44 ! Hsig 105 INTEGER, PARAMETER :: jpr_phioc = 45 ! Wave=>ocean energy flux 106 INTEGER, PARAMETER :: jpr_sdrftx = 46 ! Stokes drift on grid 1 107 INTEGER, PARAMETER :: jpr_sdrfty = 47 ! Stokes drift on grid 2 105 !** surface wave coupling ** 106 INTEGER, PARAMETER :: jpr_hsig = 44 ! Hsig 107 INTEGER, PARAMETER :: jpr_phioc = 45 ! Wave=>ocean energy flux 108 INTEGER, PARAMETER :: jpr_sdrftx = 46 ! Stokes drift on grid 1 109 INTEGER, PARAMETER :: jpr_sdrfty = 47 ! Stokes drift on grid 2 108 110 INTEGER, PARAMETER :: jpr_wper = 48 ! Mean wave period 109 111 INTEGER, PARAMETER :: jpr_wnum = 49 ! Mean wavenumber 110 INTEGER, PARAMETER :: jpr_ tauwoc= 50 ! Stress fraction adsorbed by waves112 INTEGER, PARAMETER :: jpr_wstrf = 50 ! Stress fraction adsorbed by waves 111 113 INTEGER, PARAMETER :: jpr_wdrag = 51 ! Neutral surface drag coefficient 112 INTEGER, PARAMETER :: jpr_isf = 52 113 INTEGER, PARAMETER :: jpr_icb = 53 114 INTEGER, PARAMETER :: jpr_wfreq = 54 ! Wave peak frequency 115 INTEGER, PARAMETER :: jpr_tauwx = 55 ! x component of the ocean stress from waves 116 INTEGER, PARAMETER :: jpr_tauwy = 56 ! y component of the ocean stress from waves 117 INTEGER, PARAMETER :: jpr_ts_ice = 57 ! Sea ice surface temp 118 119 INTEGER, PARAMETER :: jprcv = 57 ! total number of fields received 114 INTEGER, PARAMETER :: jpr_charn = 52 ! Chranock coefficient 115 INTEGER, PARAMETER :: jpr_twox = 53 ! wave to ocean momentum flux 116 INTEGER, PARAMETER :: jpr_twoy = 54 ! wave to ocean momentum flux 117 INTEGER, PARAMETER :: jpr_tawx = 55 ! net wave-supported stress 118 INTEGER, PARAMETER :: jpr_tawy = 56 ! net wave-supported stress 119 INTEGER, PARAMETER :: jpr_bhd = 57 ! Bernoulli head. waves' induced surface pressure 120 INTEGER, PARAMETER :: jpr_tusd = 58 ! zonal stokes transport 121 INTEGER, PARAMETER :: jpr_tvsd = 59 ! meridional stokes tranmport 122 INTEGER, PARAMETER :: jpr_isf = 60 123 INTEGER, PARAMETER :: jpr_icb = 61 124 INTEGER, PARAMETER :: jpr_ts_ice = 62 ! Sea ice surface temp 125 126 INTEGER, PARAMETER :: jprcv = 62 ! total number of fields received 120 127 121 128 INTEGER, PARAMETER :: jps_fice = 1 ! ice fraction sent to the atmosphere … … 172 179 & sn_snd_thick1, sn_snd_cond, sn_snd_mpnd , sn_snd_sstfrz, sn_snd_ttilyr 173 180 ! ! Received from the atmosphere 174 TYPE(FLD_C) :: sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_ tauw, sn_rcv_dqnsdt, sn_rcv_qsr, &181 TYPE(FLD_C) :: sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_dqnsdt, sn_rcv_qsr, & 175 182 & sn_rcv_qns , sn_rcv_emp , sn_rcv_rnf, sn_rcv_ts_ice 176 183 TYPE(FLD_C) :: sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_mslp, sn_rcv_icb, sn_rcv_isf 177 ! Send to waves184 ! ! Send to waves 178 185 TYPE(FLD_C) :: sn_snd_ifrac, sn_snd_crtw, sn_snd_wlev 179 ! Received from waves180 TYPE(FLD_C) :: sn_rcv_hsig, sn_rcv_phioc, sn_rcv_sdrfx, sn_rcv_sdrfy, sn_rcv_wper, sn_rcv_wnum, sn_rcv_tauwoc,&181 sn_rcv_wdrag, sn_rcv_wfreq186 ! ! Received from waves 187 TYPE(FLD_C) :: sn_rcv_hsig, sn_rcv_phioc, sn_rcv_sdrfx, sn_rcv_sdrfy, sn_rcv_wper, sn_rcv_wnum, & 188 & sn_rcv_wstrf, sn_rcv_wdrag, sn_rcv_charn, sn_rcv_taw, sn_rcv_bhd, sn_rcv_tusd, sn_rcv_tvsd 182 189 ! ! Other namelist parameters 183 190 INTEGER :: nn_cplmodel ! Maximum number of models to/from which NEMO is potentialy sending/receiving data … … 248 255 REAL(wp), DIMENSION(jpi,jpj) :: zacs, zaos 249 256 !! 250 NAMELIST/namsbc_cpl/ sn_snd_temp , sn_snd_alb , sn_snd_thick, sn_snd_crt , sn_snd_co2 , & 251 & sn_snd_ttilyr, sn_snd_cond , sn_snd_mpnd , sn_snd_sstfrz, sn_snd_thick1, & 252 & sn_snd_ifrac , sn_snd_crtw , sn_snd_wlev , sn_rcv_hsig , sn_rcv_phioc, & 253 & sn_rcv_w10m , sn_rcv_taumod, sn_rcv_tau , sn_rcv_dqnsdt, sn_rcv_qsr , & 254 & sn_rcv_sdrfx , sn_rcv_sdrfy , sn_rcv_wper , sn_rcv_wnum , sn_rcv_tauwoc, & 257 NAMELIST/namsbc_cpl/ sn_snd_temp , sn_snd_alb , sn_snd_thick, sn_snd_crt , sn_snd_co2 , & 258 & sn_snd_ttilyr, sn_snd_cond , sn_snd_mpnd , sn_snd_sstfrz, sn_snd_thick1, & 259 & sn_snd_ifrac , sn_snd_crtw , sn_snd_wlev , sn_rcv_hsig , sn_rcv_phioc, & 260 & sn_rcv_w10m , sn_rcv_taumod, sn_rcv_tau , sn_rcv_dqnsdt, sn_rcv_qsr , & 261 & sn_rcv_sdrfx , sn_rcv_sdrfy , sn_rcv_wper , sn_rcv_wnum , sn_rcv_wstrf, & 262 & sn_rcv_charn , sn_rcv_taw , sn_rcv_bhd , sn_rcv_tusd , sn_rcv_tvsd, & 255 263 & sn_rcv_wdrag , sn_rcv_qns , sn_rcv_emp , sn_rcv_rnf , sn_rcv_cal , & 256 264 & sn_rcv_iceflx, sn_rcv_co2 , nn_cplmodel , ln_usecplmask, sn_rcv_mslp , & 257 & sn_rcv_icb , sn_rcv_isf , sn_rcv_ wfreq , sn_rcv_tauw, nn_cats_cpl , &258 & sn_rcv_ts_ice 265 & sn_rcv_icb , sn_rcv_isf , sn_rcv_ts_ice, nn_cats_cpl 266 259 267 260 268 !!--------------------------------------------------------------------- … … 294 302 WRITE(numout,*)' sea ice heat fluxes = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')' 295 303 WRITE(numout,*)' atm co2 = ', TRIM(sn_rcv_co2%cldes ), ' (', TRIM(sn_rcv_co2%clcat ), ')' 304 WRITE(numout,*)' Sea ice surface skin temperature= ', TRIM(sn_rcv_ts_ice%cldes), ' (', TRIM(sn_rcv_ts_ice%clcat), ')' 305 WRITE(numout,*)' surface waves:' 296 306 WRITE(numout,*)' significant wave heigth = ', TRIM(sn_rcv_hsig%cldes ), ' (', TRIM(sn_rcv_hsig%clcat ), ')' 297 307 WRITE(numout,*)' wave to oce energy flux = ', TRIM(sn_rcv_phioc%cldes ), ' (', TRIM(sn_rcv_phioc%clcat ), ')' … … 300 310 WRITE(numout,*)' Mean wave period = ', TRIM(sn_rcv_wper%cldes ), ' (', TRIM(sn_rcv_wper%clcat ), ')' 301 311 WRITE(numout,*)' Mean wave number = ', TRIM(sn_rcv_wnum%cldes ), ' (', TRIM(sn_rcv_wnum%clcat ), ')' 302 WRITE(numout,*)' Wave peak frequency = ', TRIM(sn_rcv_wfreq%cldes ), ' (', TRIM(sn_rcv_wfreq%clcat ), ')' 303 WRITE(numout,*)' Stress frac adsorbed by waves = ', TRIM(sn_rcv_tauwoc%cldes), ' (', TRIM(sn_rcv_tauwoc%clcat ), ')' 304 WRITE(numout,*)' Stress components by waves = ', TRIM(sn_rcv_tauw%cldes ), ' (', TRIM(sn_rcv_tauw%clcat ), ')' 312 WRITE(numout,*)' Stress frac adsorbed by waves = ', TRIM(sn_rcv_wstrf%cldes ), ' (', TRIM(sn_rcv_wstrf%clcat ), ')' 305 313 WRITE(numout,*)' Neutral surf drag coefficient = ', TRIM(sn_rcv_wdrag%cldes ), ' (', TRIM(sn_rcv_wdrag%clcat ), ')' 306 WRITE(numout,*)' Sea ice surface skin temperature= ', TRIM(sn_rcv_ts_ice%cldes), ' (', TRIM(sn_rcv_ts_ice%clcat), ')'314 WRITE(numout,*)' Charnock coefficient = ', TRIM(sn_rcv_charn%cldes ), ' (', TRIM(sn_rcv_charn%clcat ), ')' 307 315 WRITE(numout,*)' sent fields (multiple ice categories)' 308 316 WRITE(numout,*)' surface temperature = ', TRIM(sn_snd_temp%cldes ), ' (', TRIM(sn_snd_temp%clcat ), ')' … … 605 613 cpl_wper = .TRUE. 606 614 ENDIF 607 srcv(jpr_wfreq)%clname = 'O_WFreq' ! wave peak frequency608 IF( TRIM(sn_rcv_wfreq%cldes ) == 'coupled' ) THEN609 srcv(jpr_wfreq)%laction = .TRUE.610 cpl_wfreq = .TRUE.611 ENDIF612 615 srcv(jpr_wnum)%clname = 'O_WNum' ! mean wave number 613 616 IF( TRIM(sn_rcv_wnum%cldes ) == 'coupled' ) THEN … … 615 618 cpl_wnum = .TRUE. 616 619 ENDIF 617 srcv(jpr_tauwoc)%clname = 'O_TauOce' ! stress fraction adsorbed by the wave 618 IF( TRIM(sn_rcv_tauwoc%cldes ) == 'coupled' ) THEN 619 srcv(jpr_tauwoc)%laction = .TRUE. 620 cpl_tauwoc = .TRUE. 621 ENDIF 622 srcv(jpr_tauwx)%clname = 'O_Tauwx' ! ocean stress from wave in the x direction 623 srcv(jpr_tauwy)%clname = 'O_Tauwy' ! ocean stress from wave in the y direction 624 IF( TRIM(sn_rcv_tauw%cldes ) == 'coupled' ) THEN 625 srcv(jpr_tauwx)%laction = .TRUE. 626 srcv(jpr_tauwy)%laction = .TRUE. 627 cpl_tauw = .TRUE. 620 srcv(jpr_wstrf)%clname = 'O_WStrf' ! stress fraction adsorbed by the wave 621 IF( TRIM(sn_rcv_wstrf%cldes ) == 'coupled' ) THEN 622 srcv(jpr_wstrf)%laction = .TRUE. 623 cpl_wstrf = .TRUE. 628 624 ENDIF 629 625 srcv(jpr_wdrag)%clname = 'O_WDrag' ! neutral surface drag coefficient … … 632 628 cpl_wdrag = .TRUE. 633 629 ENDIF 634 IF( srcv(jpr_tauwoc)%laction .AND. srcv(jpr_tauwx)%laction .AND. srcv(jpr_tauwy)%laction ) & 635 CALL ctl_stop( 'More than one method for modifying the ocean stress has been selected ', & 636 '(sn_rcv_tauwoc=coupled and sn_rcv_tauw=coupled)' ) 630 srcv(jpr_charn)%clname = 'O_Charn' ! neutral surface drag coefficient 631 IF( TRIM(sn_rcv_charn%cldes ) == 'coupled' ) THEN 632 srcv(jpr_charn)%laction = .TRUE. 633 cpl_charn = .TRUE. 634 ENDIF 635 srcv(jpr_bhd)%clname = 'O_Bhd' ! neutral surface drag coefficient 636 IF( TRIM(sn_rcv_bhd%cldes ) == 'coupled' ) THEN 637 srcv(jpr_bhd)%laction = .TRUE. 638 cpl_bhd = .TRUE. 639 ENDIF 640 srcv(jpr_tusd)%clname = 'O_Tusd' ! neutral surface drag coefficient 641 IF( TRIM(sn_rcv_tusd%cldes ) == 'coupled' ) THEN 642 srcv(jpr_tusd)%laction = .TRUE. 643 cpl_tusd = .TRUE. 644 ENDIF 645 srcv(jpr_tvsd)%clname = 'O_Tvsd' ! neutral surface drag coefficient 646 IF( TRIM(sn_rcv_tvsd%cldes ) == 'coupled' ) THEN 647 srcv(jpr_tvsd)%laction = .TRUE. 648 cpl_tvsd = .TRUE. 649 ENDIF 650 651 srcv(jpr_twox)%clname = 'O_Twox' ! wave to ocean momentum flux in the u direction 652 srcv(jpr_twoy)%clname = 'O_Twoy' ! wave to ocean momentum flux in the v direction 653 srcv(jpr_tawx)%clname = 'O_Tawx' ! Net wave-supported stress in the u direction 654 srcv(jpr_tawy)%clname = 'O_Tawy' ! Net wave-supported stress in the v direction 655 IF( TRIM(sn_rcv_taw%cldes ) == 'coupled' ) THEN 656 srcv(jpr_twox)%laction = .TRUE. 657 srcv(jpr_twoy)%laction = .TRUE. 658 srcv(jpr_tawx)%laction = .TRUE. 659 srcv(jpr_tawy)%laction = .TRUE. 660 cpl_taw = .TRUE. 661 ENDIF 637 662 ! 638 663 ! ! ------------------------------- ! … … 1043 1068 ENDIF 1044 1069 xcplmask(:,:,0) = 1. - SUM( xcplmask(:,:,1:nn_cplmodel), dim = 3 ) 1070 ! 1045 1071 ! 1046 1072 END SUBROUTINE sbc_cpl_init … … 1118 1144 IF( ncpl_qsr_freq /= 0) ncpl_qsr_freq = 86400 / ncpl_qsr_freq ! used by top 1119 1145 1146 IF ( ln_wave .AND. nn_components == 0 ) THEN 1147 ncpl_qsr_freq = 1; 1148 WRITE(numout,*) 'ncpl_qsr_freq is set to 1 when coupling NEMO with wave (without SAS) ' 1149 ENDIF 1120 1150 ENDIF 1121 1151 ! … … 1280 1310 IF( srcv(jpr_hsig)%laction ) hsw(:,:) = frcv(jpr_hsig)%z3(:,:,1) 1281 1311 ! 1282 ! ! ========================= !1283 ! ! Wave peak frequency !1284 ! ! ========================= !1285 IF( srcv(jpr_wfreq)%laction ) wfreq(:,:) = frcv(jpr_wfreq)%z3(:,:,1)1286 !1287 1312 ! ! ========================= ! 1288 1313 ! ! Vertical mixing Qiao ! … … 1291 1316 1292 1317 ! Calculate the 3D Stokes drift both in coupled and not fully uncoupled mode 1293 IF( srcv(jpr_sdrftx)%laction .OR. srcv(jpr_sdrfty)%laction .OR. srcv(jpr_wper)%laction&1294 .OR. srcv(jpr_hsig)%laction .OR. srcv(jpr_wfreq)%laction)THEN1318 IF( srcv(jpr_sdrftx)%laction .OR. srcv(jpr_sdrfty)%laction .OR. & 1319 srcv(jpr_wper)%laction .OR. srcv(jpr_hsig)%laction ) THEN 1295 1320 CALL sbc_stokes( Kmm ) 1296 1321 ENDIF … … 1299 1324 ! ! Stress adsorbed by waves ! 1300 1325 ! ! ========================= ! 1301 IF( srcv(jpr_tauwoc)%laction .AND. ln_tauwoc ) tauoc_wave(:,:) = frcv(jpr_tauwoc)%z3(:,:,1) 1302 1303 ! ! ========================= ! 1304 ! ! Stress component by waves ! 1305 ! ! ========================= ! 1306 IF( srcv(jpr_tauwx)%laction .AND. srcv(jpr_tauwy)%laction .AND. ln_tauw ) THEN 1307 tauw_x(:,:) = frcv(jpr_tauwx)%z3(:,:,1) 1308 tauw_y(:,:) = frcv(jpr_tauwy)%z3(:,:,1) 1309 ENDIF 1310 1326 IF( srcv(jpr_wstrf)%laction .AND. ln_tauoc ) tauoc_wave(:,:) = frcv(jpr_wstrf)%z3(:,:,1) 1327 ! 1311 1328 ! ! ========================= ! 1312 1329 ! ! Wave drag coefficient ! 1313 1330 ! ! ========================= ! 1314 1331 IF( srcv(jpr_wdrag)%laction .AND. ln_cdgw ) cdn_wave(:,:) = frcv(jpr_wdrag)%z3(:,:,1) 1315 1332 ! 1333 ! ! ========================= ! 1334 ! ! Wave drag coefficient ! 1335 ! ! ========================= ! 1336 IF( srcv(jpr_charn)%laction .AND. ln_charn ) charn(:,:) = frcv(jpr_charn)%z3(:,:,1) 1337 ! 1338 ! 1339 IF( srcv(jpr_tawx)%laction .AND. ln_taw ) tawx(:,:) = frcv(jpr_tawx)%z3(:,:,1) 1340 IF( srcv(jpr_tawy)%laction .AND. ln_taw ) tawy(:,:) = frcv(jpr_tawy)%z3(:,:,1) 1341 IF( srcv(jpr_twox)%laction .AND. ln_taw ) twox(:,:) = frcv(jpr_twox)%z3(:,:,1) 1342 IF( srcv(jpr_twoy)%laction .AND. ln_taw ) twoy(:,:) = frcv(jpr_twoy)%z3(:,:,1) 1343 ! 1344 ! ! ========================= ! 1345 ! ! wave TKE flux at sfc ! 1346 ! ! ========================= ! 1347 IF( srcv(jpr_phioc)%laction .AND. ln_phioc ) phioc(:,:) = frcv(jpr_phioc)%z3(:,:,1) 1348 ! 1349 ! ! ========================= ! 1350 ! ! Bernoulli head ! 1351 ! ! ========================= ! 1352 IF( srcv(jpr_bhd)%laction .AND. ln_bern_srfc ) bhd_wave(:,:) = frcv(jpr_bhd)%z3(:,:,1) 1353 ! 1354 ! ! ========================= ! 1355 ! ! Stokes transport u dir ! 1356 ! ! ========================= ! 1357 IF( srcv(jpr_tusd)%laction .AND. ln_breivikFV_2016 ) tusd(:,:) = frcv(jpr_tusd)%z3(:,:,1) 1358 ! 1359 ! ! ========================= ! 1360 ! ! Stokes transport v dir ! 1361 ! ! ========================= ! 1362 IF( srcv(jpr_tvsd)%laction .AND. ln_breivikFV_2016 ) tvsd(:,:) = frcv(jpr_tvsd)%z3(:,:,1) 1363 ! 1316 1364 ! Fields received by SAS when OASIS coupling 1317 1365 ! (arrays no more filled at sbcssm stage) -
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/SBC/sbcmod.F90
r12489 r12991 16 16 !! 4.0 ! 2016-06 (L. Brodeau) new general bulk formulation 17 17 !! 4.0 ! 2019-03 (F. Lemarié & G. Samson) add ABL compatibility (ln_abl=TRUE) 18 !! 4.2 ! 2020-05 (G. Madec, E. Clementi) modified wave forcing and coipling 18 19 !!---------------------------------------------------------------------- 19 20 … … 54 55 USE usrdef_sbc ! user defined: surface boundary condition 55 56 USE closea ! closed sea 57 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 56 58 ! 57 59 USE prtctl ! Print control (prt_ctl routine) … … 70 72 71 73 INTEGER :: nsbc ! type of surface boundary condition (deduced from namsbc informations) 72 74 !! * Substitutions 75 # include "do_loop_substitute.h90" 73 76 !!---------------------------------------------------------------------- 74 77 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 99 102 & nn_ice , ln_ice_embd, & 100 103 & ln_traqsr, ln_dm2dc , & 101 & ln_rnf , nn_fwb , ln_ssr , ln_apr_dyn, & 102 & ln_wave , ln_cdgw , ln_sdw , ln_tauwoc , ln_stcor , & 103 & ln_tauw , nn_lsm, nn_sdrift 104 & ln_rnf , nn_fwb , ln_ssr , ln_apr_dyn, & 105 & ln_wave , nn_lsm 104 106 !!---------------------------------------------------------------------- 105 107 ! … … 140 142 WRITE(numout,*) ' bulk formulation ln_blk = ', ln_blk 141 143 WRITE(numout,*) ' ABL formulation ln_abl = ', ln_abl 144 WRITE(numout,*) ' Surface wave (forced or coupled) ln_wave = ', ln_wave 142 145 WRITE(numout,*) ' Type of coupling (Ocean/Ice/Atmosphere) : ' 143 146 WRITE(numout,*) ' ocean-atmosphere coupled formulation ln_cpl = ', ln_cpl … … 157 160 WRITE(numout,*) ' runoff / runoff mouths ln_rnf = ', ln_rnf 158 161 WRITE(numout,*) ' nb of iterations if land-sea-mask applied nn_lsm = ', nn_lsm 159 WRITE(numout,*) ' surface wave ln_wave = ', ln_wave 160 WRITE(numout,*) ' Stokes drift corr. to vert. velocity ln_sdw = ', ln_sdw 161 WRITE(numout,*) ' vertical parametrization nn_sdrift = ', nn_sdrift 162 WRITE(numout,*) ' wave modified ocean stress ln_tauwoc = ', ln_tauwoc 163 WRITE(numout,*) ' wave modified ocean stress component ln_tauw = ', ln_tauw 164 WRITE(numout,*) ' Stokes coriolis term ln_stcor = ', ln_stcor 165 WRITE(numout,*) ' neutral drag coefficient (CORE,NCAR) ln_cdgw = ', ln_cdgw 166 ENDIF 167 ! 168 IF( .NOT.ln_wave ) THEN 169 ln_sdw = .false. ; ln_cdgw = .false. ; ln_tauwoc = .false. ; ln_tauw = .false. ; ln_stcor = .false. 170 ENDIF 171 IF( ln_sdw ) THEN 172 IF( .NOT.(nn_sdrift==jp_breivik_2014 .OR. nn_sdrift==jp_li_2017 .OR. nn_sdrift==jp_peakfr) ) & 173 CALL ctl_stop( 'The chosen nn_sdrift for Stokes drift vertical velocity must be 0, 1, or 2' ) 174 ENDIF 175 ll_st_bv2014 = ( nn_sdrift==jp_breivik_2014 ) 176 ll_st_li2017 = ( nn_sdrift==jp_li_2017 ) 177 ll_st_bv_li = ( ll_st_bv2014 .OR. ll_st_li2017 ) 178 ll_st_peakfr = ( nn_sdrift==jp_peakfr ) 179 IF( ln_tauwoc .AND. ln_tauw ) & 180 CALL ctl_stop( 'More than one method for modifying the ocean stress has been selected ', & 181 '(ln_tauwoc=.true. and ln_tauw=.true.)' ) 182 IF( ln_tauwoc ) & 183 CALL ctl_warn( 'You are subtracting the wave stress to the ocean (ln_tauwoc=.true.)' ) 184 IF( ln_tauw ) & 185 CALL ctl_warn( 'The wave modified ocean stress components are used (ln_tauw=.true.) ', & 186 'This will override any other specification of the ocean stress' ) 162 ENDIF 187 163 ! 188 164 IF( .NOT.ln_usr ) THEN ! the model calendar needs some specificities (except in user defined case) … … 364 340 IF( nn_ice == 3 ) CALL cice_sbc_init( nsbc, Kbb, Kmm ) ! CICE initialization 365 341 ! 366 IF( ln_wave ) CALL sbc_wave_init ! surface wave initialisation 342 IF( ln_wave ) THEN 343 CALL sbc_wave_init ! surface wave initialisation 344 ELSE 345 IF(lwp) WRITE(numout,*) 346 IF(lwp) WRITE(numout,*) ' No surface waves : all wave related logical set to false' 347 ln_sdw = .false. 348 ln_stcor = .false. 349 ln_cdgw = .false. 350 ln_tauoc = .false. 351 ln_wave_test = .false. 352 ln_charn = .false. 353 ln_taw = .false. 354 ln_phioc = .false. 355 ln_bern_srfc = .false. 356 ln_breivikFV_2016 = .false. 357 ln_vortex_force = .false. 358 ln_stshear = .false. 359 ENDIF 367 360 ! 368 361 IF( lwxios ) THEN … … 397 390 INTEGER, INTENT(in) :: kt ! ocean time step 398 391 INTEGER, INTENT(in) :: Kbb, Kmm ! ocean time level indices 392 INTEGER :: jj, ji ! dummy loop argument 399 393 ! 400 394 LOGICAL :: ll_sas, ll_opa ! local logical … … 429 423 ! 430 424 IF( .NOT.ll_sas ) CALL sbc_ssm ( kt, Kbb, Kmm ) ! mean ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m) 431 IF( ln_wave ) CALL sbc_wave( kt, Kmm ) ! surface waves432 433 425 ! 434 426 ! !== sbc formulation ==! 435 427 ! 428 ! 436 429 SELECT CASE( nsbc ) ! Compute ocean surface boundary condition 437 430 ! ! (i.e. utau,vtau, qns, qsr, emp, sfx) … … 440 433 CASE( jp_blk ) 441 434 IF( ll_sas ) CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice, Kbb, Kmm ) ! OPA-SAS coupling: SAS receiving fields from OPA 435 !!!!!!!!!!! ATTENTION:ln_wave is not only used for oasis coupling !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 436 IF( ln_wave ) THEN 437 IF ( lk_oasis ) CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice, Kbb, Kmm ) ! OPA-wave coupling 438 CALL sbc_wave ( kt, Kmm ) 439 ENDIF 442 440 CALL sbc_blk ( kt ) ! bulk formulation for the ocean 443 441 ! … … 453 451 IF( ln_mixcpl ) CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice, Kbb, Kmm ) ! forced-coupled mixed formulation after forcing 454 452 ! 455 IF ( ln_wave .AND. (ln_tauwoc .OR. ln_tauw) ) CALL sbc_wstress( ) ! Wind stress provided by waves 453 IF( ln_wave .AND. ln_tauoc ) THEN ! Wave stress reduction 454 DO_2D_00_00 455 utau(ji,jj) = utau(ji,jj) * ( tauoc_wave(ji,jj) + tauoc_wave(ji-1,jj) ) * 0.5_wp 456 vtau(ji,jj) = vtau(ji,jj) * ( tauoc_wave(ji,jj) + tauoc_wave(ji,jj-1) ) * 0.5_wp 457 END_2D 458 ! 459 CALL lbc_lnk( 'sbcwave', utau, 'U', -1. ) 460 CALL lbc_lnk( 'sbcwave', vtau, 'V', -1. ) 461 ! 462 taum(:,:) = taum(:,:)*tauoc_wave(:,:) 463 ! 464 IF( kt == nit000 ) CALL ctl_warn( 'sbc: You are subtracting the wave stress to the ocean.', & 465 & 'If not requested select ln_tauoc=.false.' ) 466 ! 467 ELSEIF( ln_wave .AND. ln_taw ) THEN ! Wave stress reduction 468 utau(:,:) = utau(:,:) - tawx(:,:) + twox(:,:) 469 vtau(:,:) = vtau(:,:) - tawy(:,:) + twoy(:,:) 470 CALL lbc_lnk( 'sbcwave', utau, 'U', -1. ) 471 CALL lbc_lnk( 'sbcwave', vtau, 'V', -1. ) 472 ! 473 DO_2D_00_00 474 taum(ji,jj) = sqrt((.5*(utau(ji-1,jj)+utau(ji,jj)))**2 + (.5*(vtau(ji,jj-1)+vtau(ji,jj)))**2) 475 END_2D 476 ! 477 IF( kt == nit000 ) CALL ctl_warn( 'sbc: You are subtracting the wave stress to the ocean.', & 478 & 'If not requested select ln_taw=.false.' ) 479 ! 480 ENDIF 481 CALL lbc_lnk( 'sbcmod', taum(:,:), 'T', 1. ) 456 482 ! 457 483 ! !== Misc. Options ==! -
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/SBC/sbcwave.F90
r12377 r12991 9 9 !! - ! 2016-12 (G. Madec, E. Clementi) update Stoke drift computation 10 10 !! + add sbc_wave_ini routine 11 !! 4.2 ! 2020-06 (G. Madec, E. Clement) updates, new Stoke drift computation 12 !! according to Covelard et al.,2019 11 13 !!---------------------------------------------------------------------- 12 14 13 15 !!---------------------------------------------------------------------- 14 16 !! sbc_stokes : calculate 3D Stokes-drift velocities 15 !! sbc_wave : wave data from wave model in netcdf files17 !! sbc_wave : wave data from wave model: forced (netcdf files) or coupled mode 16 18 !! sbc_wave_init : initialisation fo surface waves 17 19 !!---------------------------------------------------------------------- 18 USE phycst ! physical constants 20 USE phycst ! physical constants 19 21 USE oce ! ocean variables 20 USE sbc_oce ! Surface boundary condition: ocean fields21 USE zdf_oce, ONLY : ln_zdfswm22 USE dom_oce ! ocean domain variables 23 USE sbc_oce ! Surface boundary condition: ocean fields 22 24 USE bdy_oce ! open boundary condition variables 23 25 USE domvvl ! domain: variable volume layers … … 26 28 USE in_out_manager ! I/O manager 27 29 USE lib_mpp ! distribued memory computing library 28 USE fldread 30 USE fldread ! read input fields 29 31 30 32 IMPLICIT NONE … … 32 34 33 35 PUBLIC sbc_stokes ! routine called in sbccpl 34 PUBLIC sbc_wstress ! routine called in sbcmod35 36 PUBLIC sbc_wave ! routine called in sbcmod 36 37 PUBLIC sbc_wave_init ! routine called in sbcmod 37 38 38 39 ! Variables checking if the wave parameters are coupled (if not, they are read from file) 39 LOGICAL, PUBLIC :: cpl_hsig = .FALSE. 40 LOGICAL, PUBLIC :: cpl_phioc = .FALSE. 41 LOGICAL, PUBLIC :: cpl_sdrftx = .FALSE. 42 LOGICAL, PUBLIC :: cpl_sdrfty = .FALSE. 43 LOGICAL, PUBLIC :: cpl_wper = .FALSE. 44 LOGICAL, PUBLIC :: cpl_wfreq = .FALSE. 45 LOGICAL, PUBLIC :: cpl_wnum = .FALSE. 46 LOGICAL, PUBLIC :: cpl_tauwoc = .FALSE. 47 LOGICAL, PUBLIC :: cpl_tauw = .FALSE. 48 LOGICAL, PUBLIC :: cpl_wdrag = .FALSE. 40 LOGICAL, PUBLIC :: cpl_hsig = .FALSE. 41 LOGICAL, PUBLIC :: cpl_phioc = .FALSE. 42 LOGICAL, PUBLIC :: cpl_sdrftx = .FALSE. 43 LOGICAL, PUBLIC :: cpl_sdrfty = .FALSE. 44 LOGICAL, PUBLIC :: cpl_wper = .FALSE. 45 LOGICAL, PUBLIC :: cpl_wnum = .FALSE. 46 LOGICAL, PUBLIC :: cpl_wstrf = .FALSE. 47 LOGICAL, PUBLIC :: cpl_wdrag = .FALSE. 48 LOGICAL, PUBLIC :: cpl_charn = .FALSE. 49 LOGICAL, PUBLIC :: cpl_taw = .FALSE. 50 LOGICAL, PUBLIC :: cpl_bhd = .FALSE. 51 LOGICAL, PUBLIC :: cpl_tusd = .FALSE. 52 LOGICAL, PUBLIC :: cpl_tvsd = .FALSE. 49 53 50 54 INTEGER :: jpfld ! number of files to read for stokes drift … … 53 57 INTEGER :: jp_hsw ! index of significant wave hight (m) at T-point 54 58 INTEGER :: jp_wmp ! index of mean wave period (s) at T-point 55 INTEGER :: jp_wfr ! index of wave peak frequency (1/s) at T-point56 59 57 60 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_cd ! structure of input fields (file informations, fields read) Drag Coefficient 58 61 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_sd ! structure of input fields (file informations, fields read) Stokes Drift 59 62 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_wn ! structure of input fields (file informations, fields read) wave number for Qiao 60 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_tauw ! structure of input fields (file informations, fields read) ocean stress components from wave model 62 63 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: cdn_wave !: 64 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: hsw, wmp, wnum !: 65 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: wfreq !: 66 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: tauoc_wave !: 67 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: tauw_x, tauw_y !: 68 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: tsd2d !: 69 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: div_sd !: barotropic stokes drift divergence 70 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: ut0sd, vt0sd !: surface Stokes drift velocities at t-point 71 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: usd , vsd , wsd !: Stokes drift velocities at u-, v- & w-points, resp. 72 63 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_tauoc ! structure of input fields (file informations, fields read) normalized wave stress into the ocean 64 65 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: cdn_wave !: Neutral drag coefficient at t-point 66 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: hsw !: Significant Wave Height at t-point 67 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: wmp !: Wave Mean Period at t-point 68 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: wnum !: Wave Number at t-point 69 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: tauoc_wave !: stress reduction factor at t-point 70 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: tsd2d !: Surface Stokes Drift module at t-point 71 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: div_sd !: barotropic stokes drift divergence 72 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: ut0sd, vt0sd !: surface Stokes drift velocities at t-point 73 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: usd, vsd, wsd !: Stokes drift velocities at u-, v- & w-points, resp.u 74 ! 75 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: charn !: charnock coefficient at t-point 76 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: tawx !: Net wave-supported stress, u 77 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: tawy !: Net wave-supported stress, v 78 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: twox !: wave-ocean momentum flux, u 79 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: twoy !: wave-ocean momentum flux, v 80 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: tauoc_wavex !: stress reduction factor at, u component 81 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: tauoc_wavey !: stress reduction factor at, v component 82 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: phioc !: tke flux from wave model 83 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: KZN2 !: Kz*N2 84 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: bhd_wave !: Bernoulli head. wave induce pression 85 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: tusd, tvsd !: Stokes drift transport 86 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: ZMX !: Kz*N2 73 87 !! * Substitutions 74 88 # include "do_loop_substitute.h90" … … 87 101 !! 2014 (DOI: 10.1175/JPO-D-14-0020.1) 88 102 !! 89 !! ** Method : - Calculate Stokes transport speed 90 !! - Calculate horizontal divergence 91 !! - Integrate the horizontal divergenze from the bottom 92 !! ** action 103 !! ** Method : - Calculate the horizontal Stokes drift velocity (Breivik et al. 2014) 104 !! - Calculate its horizontal divergence 105 !! - Calculate the vertical Stokes drift velocity 106 !! - Calculate the barotropic Stokes drift divergence 107 !! 108 !! ** action : - tsd2d : module of the surface Stokes drift velocity 109 !! - usd, vsd, wsd : 3 components of the Stokes drift velocity 110 !! - div_sd : barotropic Stokes drift divergence 93 111 !!--------------------------------------------------------------------- 94 112 INTEGER, INTENT(in) :: Kmm ! ocean time level index 95 113 INTEGER :: jj, ji, jk ! dummy loop argument 96 114 INTEGER :: ik ! local integer 97 REAL(wp) :: ztransp, zfac, zsp0 98 REAL(wp) :: zdepth, zsqrt_depth, zexp_depth, z_two_thirds, zsqrtpi !sqrt of pi 99 REAL(wp) :: zbot_u, zbot_v, zkb_u, zkb_v, zke3_u, zke3_v, zda_u, zda_v 100 REAL(wp) :: zstokes_psi_u_bot, zstokes_psi_v_bot 101 REAL(wp) :: zdep_u, zdep_v, zkh_u, zkh_v 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) ) 115 REAL(wp) :: ztransp, zfac, ztemp, zsp0, zsqrt, zbreiv16_w 116 REAL(wp) :: zdep_u, zdep_v, zkh_u, zkh_v, zda_u, zda_v, sdtrp 117 REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: zk_t, zk_u, zk_v, zu0_sd, zv0_sd ! 2D workspace 118 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ze3divh, zInt_w ! 3D workspace 119 !!--------------------------------------------------------------------- 120 ! 121 ALLOCATE( ze3divh(jpi,jpj,jpk), zInt_w(jpi,jpj,jpk)) 108 122 ALLOCATE( zk_t(jpi,jpj), zk_u(jpi,jpj), zk_v(jpi,jpj), zu0_sd(jpi,jpj), zv0_sd(jpi,jpj) ) 123 zk_t (:,:) = 0._wp 124 zk_u (:,:) = 0._wp 125 zk_v (:,:) = 0._wp 126 zu0_sd (:,:) = 0._wp 127 zv0_sd (:,:) = 0._wp 128 ze3divh (:,:,:) = 0._wp 129 109 130 ! 110 131 ! select parameterization for the calculation of vertical Stokes drift 111 132 ! exp. wave number at t-point 112 IF( ll_st_bv_li ) THEN ! (Eq. (19) in Breivik et al. (2014) ) 133 IF( ln_breivikFV_2016 ) THEN 134 ! Assumptions : ut0sd and vt0sd are surface Stokes drift at T-points 135 ! sdtrp is the norm of Stokes transport 136 ! 137 zfac = 0.166666666667_wp 138 DO_2D_11_11 ! In the deep-water limit we have ke = ||ust0||/( 6 * ||transport|| ) 139 zsp0 = SQRT( ut0sd(ji,jj)*ut0sd(ji,jj) + vt0sd(ji,jj)*vt0sd(ji,jj) ) !<-- norm of Surface Stokes drift 140 tsd2d(ji,jj) = zsp0 141 IF( cpl_tusd .AND. cpl_tvsd ) THEN !stokes transport is provided in coupled mode 142 sdtrp = SQRT( tusd(ji,jj)*tusd(ji,jj) + tvsd(ji,jj)*tvsd(ji,jj) ) !<-- norm of Surface Stokes drift transport 143 ELSE 144 ! Stokes drift transport estimated from Hs and Tmean 145 sdtrp = 2.0_wp * rpi / 16.0_wp * & 146 & hsw(ji,jj)*hsw(ji,jj) / MAX( wmp(ji,jj), 0.0000001_wp ) 147 ENDIF 148 zk_t (ji,jj) = zfac * zsp0 / MAX ( sdtrp, 0.0000001_wp ) !<-- ke = ||ust0||/( 6 * ||transport|| ) 149 END_2D 150 !# define zInt_w ze3divh 151 DO_3D_11_11( 1, jpk ) ! Compute the primitive of Breivik 2016 function at W-points 152 zfac = - 2._wp * zk_t (ji,jj) * gdepw(ji,jj,jk,Kmm) !<-- zfac should be negative definite 153 ztemp = EXP ( zfac ) 154 zsqrt = SQRT( -zfac ) 155 zbreiv16_w = ztemp - SQRT(rpi)*zsqrt*ERFC(zsqrt) !Eq. 16 Breivik 2016 156 zInt_w(ji,jj,jk) = ztemp - 4._wp * zk_t (ji,jj) * gdepw(ji,jj,jk,Kmm) * zbreiv16_w 157 END_3D 158 ! 159 DO jk = 1, jpkm1 160 zfac = 0.166666666667_wp 161 DO_2D_11_11 !++ Compute the FV Breivik 2016 function at T-points 162 zsp0 = zfac / MAX(zk_t (ji,jj),0.0000001_wp) 163 ztemp = zInt_w(ji,jj,jk) - zInt_w(ji,jj,jk+1) 164 zu0_sd(ji,jj) = ut0sd(ji,jj) * zsp0 * ztemp * tmask(ji,jj,jk) 165 zv0_sd(ji,jj) = vt0sd(ji,jj) * zsp0 * ztemp * tmask(ji,jj,jk) 166 END_2D 167 DO_2D_10_10 ! ++ Interpolate at U/V points 168 zfac = 1.0_wp / e3u(ji ,jj,jk,Kmm) 169 usd(ji,jj,jk) = 0.5_wp * zfac * ( zu0_sd(ji,jj)+zu0_sd(ji+1,jj) ) * umask(ji,jj,jk) 170 zfac = 1.0_wp / e3v(ji ,jj,jk,Kmm) 171 vsd(ji,jj,jk) = 0.5_wp * zfac * ( zv0_sd(ji,jj)+zv0_sd(ji,jj+1) ) * vmask(ji,jj,jk) 172 END_2D 173 ENDDO 174 !# undef zInt_w 175 ! 176 ELSE 113 177 zfac = 2.0_wp * rpi / 16.0_wp 114 178 DO_2D_11_11 … … 127 191 zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) ) 128 192 END_2D 129 ELSE IF( ll_st_peakfr ) THEN ! peak wave number calculated from the peak frequency received by the wave model 130 DO_2D_11_11 131 zk_t(ji,jj) = ( 2.0_wp * rpi * wfreq(ji,jj) ) * ( 2.0_wp * rpi * wfreq(ji,jj) ) / grav 132 END_2D 133 DO_2D_10_10 134 zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) ) 135 zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) ) 136 ! 137 zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) ) 138 zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) ) 139 END_2D 140 ENDIF 141 ! 193 142 194 ! !== horizontal Stokes Drift 3D velocity ==! 143 IF( ll_st_bv2014 ) THEN 195 144 196 DO_3D_00_00( 1, jpkm1 ) 145 197 zdep_u = 0.5_wp * ( gdept(ji,jj,jk,Kmm) + gdept(ji+1,jj,jk,Kmm) ) 146 198 zdep_v = 0.5_wp * ( gdept(ji,jj,jk,Kmm) + gdept(ji,jj+1,jk,Kmm) ) 147 ! 199 ! 148 200 zkh_u = zk_u(ji,jj) * zdep_u ! k * depth 149 201 zkh_v = zk_v(ji,jj) * zdep_v … … 155 207 vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk) 156 208 END_3D 157 ELSE IF( ll_st_li2017 .OR. ll_st_peakfr ) THEN158 ALLOCATE( zstokes_psi_u_top(jpi,jpj), zstokes_psi_v_top(jpi,jpj) )159 DO_2D_10_10160 zstokes_psi_u_top(ji,jj) = 0._wp161 zstokes_psi_v_top(ji,jj) = 0._wp162 END_2D163 zsqrtpi = SQRT(rpi)164 z_two_thirds = 2.0_wp / 3.0_wp165 DO_3D_00_00( 1, jpkm1 )166 zbot_u = ( gdepw(ji,jj,jk+1,Kmm) + gdepw(ji+1,jj,jk+1,Kmm) ) ! 2 * bottom depth167 zbot_v = ( gdepw(ji,jj,jk+1,Kmm) + gdepw(ji,jj+1,jk+1,Kmm) ) ! 2 * bottom depth168 zkb_u = zk_u(ji,jj) * zbot_u ! 2 * k * bottom depth169 zkb_v = zk_v(ji,jj) * zbot_v ! 2 * k * bottom depth170 !171 zke3_u = MAX(1.e-8_wp, 2.0_wp * zk_u(ji,jj) * e3u(ji,jj,jk,Kmm)) ! 2k * thickness172 zke3_v = MAX(1.e-8_wp, 2.0_wp * zk_v(ji,jj) * e3v(ji,jj,jk,Kmm)) ! 2k * thickness173 174 ! Depth attenuation .... do u component first..175 zdepth = zkb_u176 zsqrt_depth = SQRT(zdepth)177 zexp_depth = EXP(-zdepth)178 zstokes_psi_u_bot = 1.0_wp - zexp_depth &179 & - z_two_thirds * ( zsqrtpi*zsqrt_depth*zdepth*ERFC(zsqrt_depth) &180 & + 1.0_wp - (1.0_wp + zdepth)*zexp_depth )181 zda_u = ( zstokes_psi_u_bot - zstokes_psi_u_top(ji,jj) ) / zke3_u182 zstokes_psi_u_top(ji,jj) = zstokes_psi_u_bot183 184 ! ... and then v component185 zdepth =zkb_v186 zsqrt_depth = SQRT(zdepth)187 zexp_depth = EXP(-zdepth)188 zstokes_psi_v_bot = 1.0_wp - zexp_depth &189 & - z_two_thirds * ( zsqrtpi*zsqrt_depth*zdepth*ERFC(zsqrt_depth) &190 & + 1.0_wp - (1.0_wp + zdepth)*zexp_depth )191 zda_v = ( zstokes_psi_v_bot - zstokes_psi_v_top(ji,jj) ) / zke3_v192 zstokes_psi_v_top(ji,jj) = zstokes_psi_v_bot193 !194 usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk)195 vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk)196 END_3D197 DEALLOCATE( zstokes_psi_u_top, zstokes_psi_v_top )198 209 ENDIF 199 210 … … 242 253 CALL iom_put( "vstokes", vsd ) 243 254 CALL iom_put( "wstokes", wsd ) 244 !245 DEALLOCATE( ze3divh )255 ! ! 256 DEALLOCATE( ze3divh, zInt_w ) 246 257 DEALLOCATE( zk_t, zk_u, zk_v, zu0_sd, zv0_sd ) 247 258 ! 248 259 END SUBROUTINE sbc_stokes 249 250 251 SUBROUTINE sbc_wstress( ) 252 !!--------------------------------------------------------------------- 253 !! *** ROUTINE sbc_wstress *** 254 !! 255 !! ** Purpose : Updates the ocean momentum modified by waves 256 !! 257 !! ** Method : - Calculate u,v components of stress depending on stress 258 !! model 259 !! - Calculate the stress module 260 !! - The wind module is not modified by waves 261 !! ** action 262 !!--------------------------------------------------------------------- 263 INTEGER :: jj, ji ! dummy loop argument 264 ! 265 IF( ln_tauwoc ) THEN 266 utau(:,:) = utau(:,:)*tauoc_wave(:,:) 267 vtau(:,:) = vtau(:,:)*tauoc_wave(:,:) 268 taum(:,:) = taum(:,:)*tauoc_wave(:,:) 269 ENDIF 270 ! 271 IF( ln_tauw ) THEN 272 DO_2D_10_10 273 ! Stress components at u- & v-points 274 utau(ji,jj) = 0.5_wp * ( tauw_x(ji,jj) + tauw_x(ji+1,jj) ) 275 vtau(ji,jj) = 0.5_wp * ( tauw_y(ji,jj) + tauw_y(ji,jj+1) ) 276 ! 277 ! Stress module at t points 278 taum(ji,jj) = SQRT( tauw_x(ji,jj)*tauw_x(ji,jj) + tauw_y(ji,jj)*tauw_y(ji,jj) ) 279 END_2D 280 CALL lbc_lnk_multi( 'sbcwave', utau(:,:), 'U', -1. , vtau(:,:), 'V', -1. , taum(:,:) , 'T', -1. ) 281 ENDIF 282 ! 283 END SUBROUTINE sbc_wstress 284 285 260 ! 261 ! 286 262 SUBROUTINE sbc_wave( kt, Kmm ) 287 263 !!--------------------------------------------------------------------- 288 264 !! *** ROUTINE sbc_wave *** 289 265 !! 290 !! ** Purpose : read wave parameters from wave model in netcdf files. 291 !! 292 !! ** Method : - Read namelist namsbc_wave 293 !! - Read Cd_n10 fields in netcdf files 294 !! - Read stokes drift 2d in netcdf files 295 !! - Read wave number in netcdf files 296 !! - Compute 3d stokes drift using Breivik et al.,2014 297 !! formulation 298 !! ** action 266 !! ** Purpose : read wave parameters from wave model in netcdf files 267 !! or from a coupled wave mdoel 268 !! 299 269 !!--------------------------------------------------------------------- 300 270 INTEGER, INTENT(in ) :: kt ! ocean time step 301 271 INTEGER, INTENT(in ) :: Kmm ! ocean time index 302 272 !!--------------------------------------------------------------------- 273 ! 274 IF( kt == nit000 .AND. lwp ) THEN 275 WRITE(numout,*) 276 WRITE(numout,*) 'sbc_wave : update the read waves fields' 277 WRITE(numout,*) '~~~~~~~~ ' 278 ENDIF 303 279 ! 304 280 IF( ln_cdgw .AND. .NOT. cpl_wdrag ) THEN !== Neutral drag coefficient ==! … … 307 283 ENDIF 308 284 309 IF( ln_tauwoc .AND. .NOT. cpl_tauwoc ) THEN !== Wave induced stress ==! 310 CALL fld_read( kt, nn_fsbc, sf_tauwoc ) ! read wave norm stress from external forcing 311 tauoc_wave(:,:) = sf_tauwoc(1)%fnow(:,:,1) * tmask(:,:,1) 312 ENDIF 313 314 IF( ln_tauw .AND. .NOT. cpl_tauw ) THEN !== Wave induced stress ==! 315 CALL fld_read( kt, nn_fsbc, sf_tauw ) ! read ocean stress components from external forcing (T grid) 316 tauw_x(:,:) = sf_tauw(1)%fnow(:,:,1) * tmask(:,:,1) 317 tauw_y(:,:) = sf_tauw(2)%fnow(:,:,1) * tmask(:,:,1) 318 ENDIF 319 320 IF( ln_sdw ) THEN !== Computation of the 3d Stokes Drift ==! 285 IF( ln_tauoc .AND. .NOT. cpl_wstrf ) THEN !== Wave induced stress ==! 286 CALL fld_read( kt, nn_fsbc, sf_tauoc ) ! read stress reduction factor due to wave from external forcing 287 tauoc_wave(:,:) = sf_tauoc(1)%fnow(:,:,1) * tmask(:,:,1) 288 ELSEIF ( ln_taw .AND. cpl_taw ) THEN 289 IF (kt < 1) THEN ! The first fields gave by OASIS have very high erroneous values .... 290 twox(:,:)=0._wp 291 twoy(:,:)=0._wp 292 tawx(:,:)=0._wp 293 tawy(:,:)=0._wp 294 tauoc_wavex(:,:) = 1._wp 295 tauoc_wavey(:,:) = 1._wp 296 ELSE 297 tauoc_wavex(:,:) = abs(twox(:,:)/tawx(:,:)) 298 tauoc_wavey(:,:) = abs(twoy(:,:)/tawy(:,:)) 299 ENDIF 300 ENDIF 301 302 IF ( ln_phioc .and. cpl_phioc .and. kt == nit000 ) THEN 303 WRITE(numout,*) 304 WRITE(numout,*) 'sbc_wave : PHIOC from wave model' 305 WRITE(numout,*) '~~~~~~~~ ' 306 ENDIF 307 308 IF( ln_sdw .AND. .NOT. cpl_sdrftx) THEN !== Computation of the 3d Stokes Drift ==! 321 309 ! 322 310 IF( jpfld > 0 ) THEN ! Read from file only if the field is not coupled 323 311 CALL fld_read( kt, nn_fsbc, sf_sd ) ! read wave parameters from external forcing 312 ! ! NB: test case mode, not read as jpfld=0 324 313 IF( jp_hsw > 0 ) hsw (:,:) = sf_sd(jp_hsw)%fnow(:,:,1) * tmask(:,:,1) ! significant wave height 325 314 IF( jp_wmp > 0 ) wmp (:,:) = sf_sd(jp_wmp)%fnow(:,:,1) * tmask(:,:,1) ! wave mean period 326 IF( jp_wfr > 0 ) wfreq(:,:) = sf_sd(jp_wfr)%fnow(:,:,1) * tmask(:,:,1) ! Peak wave frequency327 315 IF( jp_usd > 0 ) ut0sd(:,:) = sf_sd(jp_usd)%fnow(:,:,1) * tmask(:,:,1) ! 2D zonal Stokes Drift at T point 328 316 IF( jp_vsd > 0 ) vt0sd(:,:) = sf_sd(jp_vsd)%fnow(:,:,1) * tmask(:,:,1) ! 2D meridional Stokes Drift at T point 329 317 ENDIF 330 318 ! 331 ! Read also wave number if needed, so that it is available in coupling routines 332 IF( ln_zdfswm .AND. .NOT.cpl_wnum ) THEN 333 CALL fld_read( kt, nn_fsbc, sf_wn ) ! read wave parameters from external forcing 334 wnum(:,:) = sf_wn(1)%fnow(:,:,1) * tmask(:,:,1) 335 ENDIF 336 337 ! Calculate only if required fields have been read 338 ! In coupled wave model-NEMO case the call is done after coupling 319 IF( jpfld == 4 .OR. ln_wave_test ) & 320 & CALL sbc_stokes( kt ) ! Calculate only if all required fields are read 321 ! ! or in wave test case 322 ! ! ! In coupled case the call is done after (in sbc_cpl) 323 ENDIF 339 324 ! 340 IF( ( ll_st_bv_li .AND. jp_hsw>0 .AND. jp_wmp>0 .AND. jp_usd>0 .AND. jp_vsd>0 ) .OR. &341 & ( ll_st_peakfr .AND. jp_wfr>0 .AND. jp_usd>0 .AND. jp_vsd>0 ) ) CALL sbc_stokes( Kmm )342 !343 ENDIF344 !345 325 END SUBROUTINE sbc_wave 346 326 … … 350 330 !! *** ROUTINE sbc_wave_init *** 351 331 !! 352 !! ** Purpose : read wave parameters from wave model in netcdf files.332 !! ** Purpose : Initialisation fo surface waves 353 333 !! 354 334 !! ** Method : - Read namelist namsbc_wave 355 !! - Read Cd_n10 fields in netcdf files 356 !! - Read stokes drift 2d in netcdf files 357 !! - Read wave number in netcdf files 358 !! - Compute 3d stokes drift using Breivik et al.,2014 359 !! formulation 335 !! - create the structure used to read required wave fields 336 !! (its size depends on namelist options) 360 337 !! ** action 361 338 !!--------------------------------------------------------------------- … … 364 341 !! 365 342 CHARACTER(len=100) :: cn_dir ! Root directory for location of drag coefficient files 366 TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) :: slf_i , slf_j! array of namelist informations on the fields to read343 TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) :: slf_i ! array of namelist informations on the fields to read 367 344 TYPE(FLD_N) :: sn_cdg, sn_usd, sn_vsd, & 368 & sn_hsw, sn_wmp, sn_wfr, sn_wnum, & 369 & sn_tauwoc, sn_tauwx, sn_tauwy ! informations about the fields to be read 370 ! 371 NAMELIST/namsbc_wave/ sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wfr, & 372 sn_wnum, sn_tauwoc, sn_tauwx, sn_tauwy 373 !!--------------------------------------------------------------------- 345 & sn_hsw, sn_wmp, sn_wnum, sn_tauoc ! informations about the fields to be read 346 ! 347 NAMELIST/namsbc_wave/ cn_dir, sn_cdg, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wnum, sn_tauoc, & 348 & ln_cdgw, ln_sdw, ln_tauoc, ln_stcor, ln_charn, ln_taw, ln_phioc, & 349 & ln_wave_test, ln_bern_srfc, ln_breivikFV_2016, ln_vortex_force, ln_stshear 350 !!--------------------------------------------------------------------- 351 IF(lwp) THEN 352 WRITE(numout,*) 353 WRITE(numout,*) 'sbc_wave_init : surface waves in the system' 354 WRITE(numout,*) '~~~~~~~~~~~~~ ' 355 ENDIF 374 356 ! 375 357 READ ( numnam_ref, namsbc_wave, IOSTAT = ios, ERR = 901) 376 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in reference namelist' 377 358 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in reference namelist') 359 378 360 READ ( numnam_cfg, namsbc_wave, IOSTAT = ios, ERR = 902 ) 379 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namsbc_wave in configuration 361 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namsbc_wave in configurationnamelist' ) 380 362 IF(lwm) WRITE ( numond, namsbc_wave ) 381 363 ! 382 IF( ln_cdgw ) THEN 383 IF( .NOT. cpl_wdrag ) THEN 384 ALLOCATE( sf_cd(1), STAT=ierror ) !* allocate and fill sf_wave with sn_cdg 385 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' ) 364 IF(lwp) THEN 365 WRITE(numout,*) ' Namelist namsbc_wave' 366 WRITE(numout,*) ' Stokes drift ln_sdw = ', ln_sdw 367 WRITE(numout,*) ' Stokes Coriolis & tracer advection terms ln_stcor = ', ln_stcor 368 WRITE(numout,*) ' wave modified ocean stress ln_tauoc = ', ln_tauoc 369 WRITE(numout,*) ' neutral drag coefficient (CORE bulk only) ln_cdgw = ', ln_cdgw 370 WRITE(numout,*) ' charnock coefficient (IFS bulk only) ln_charn = ', ln_charn 371 WRITE(numout,*) ' Test with constant wave fields ln_wave_test = ', ln_wave_test 372 ENDIF 373 374 ! ! option check 375 IF( .NOT.( ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor .OR. ln_charn) ) & 376 & CALL ctl_warn( 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauoc=F, ln_stcor=F') 377 IF( ln_cdgw .AND. ln_blk ) & 378 & CALL ctl_stop( 'drag coefficient read from wave model NOT available yet with aerobulk package') 379 IF( ln_stcor .AND. .NOT.ln_sdw ) & 380 & CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T') 381 382 ! !== Allocate wave arrays ==! 383 ALLOCATE( ut0sd (jpi,jpj) , vt0sd (jpi,jpj) ) 384 ALLOCATE( hsw (jpi,jpj) , wmp (jpi,jpj) ) 385 ALLOCATE( wnum (jpi,jpj) ) 386 ALLOCATE( tsd2d (jpi,jpj) , div_sd(jpi,jpj) , bhd_wave(jpi,jpj) ) 387 ALLOCATE( usd (jpi,jpj,jpk), vsd (jpi,jpj,jpk), wsd (jpi,jpj,jpk) ) 388 ALLOCATE( tusd (jpi,jpj) , tvsd (jpi,jpj) , ZMX (jpi,jpj,jpk) ) 389 usd (:,:,:) = 0._wp 390 vsd (:,:,:) = 0._wp 391 wsd (:,:,:) = 0._wp 392 hsw (:,:) = 0._wp 393 wmp (:,:) = 0._wp 394 ut0sd (:,:) = 0._wp 395 vt0sd (:,:) = 0._wp 396 tusd (:,:) = 0._wp 397 tvsd (:,:) = 0._wp 398 bhd_wave(:,:) = 0._wp 399 ZMX (:,:,:) = 0._wp 400 ! 401 IF( ln_wave_test ) THEN !== Wave TEST case ==! set uniform waves fields 402 jpfld = 0 ! No field read 403 ln_cdgw = .FALSE. ! No neutral wave drag input 404 ln_tauoc = .FALSE. ! No wave induced drag reduction factor 405 ut0sd(:,:) = 0.13_wp * tmask(:,:,1) ! m/s 406 vt0sd(:,:) = 0.00_wp ! m/s 407 hsw (:,:) = 2.80_wp ! meters 408 wmp (:,:) = 8.00_wp ! seconds 409 ! 410 ELSE !== create the structure associated with fields to be read ==! 411 IF( ln_cdgw ) THEN ! wave drag 412 IF( .NOT. cpl_wdrag ) THEN 413 ALLOCATE( sf_cd(1), STAT=ierror ) !* allocate and fill sf_wave with sn_cdg 414 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' ) 415 ! 416 ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1) ) 417 IF( sn_cdg%ln_tint ) ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) ) 418 CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave_init', 'Wave module ', 'namsbc_wave' ) 419 ENDIF 420 ALLOCATE( cdn_wave(jpi,jpj) ) 421 cdn_wave(:,:) = 0._wp 422 ENDIF 423 IF( ln_charn ) THEN ! wave drag 424 IF( .NOT. cpl_charn ) THEN 425 CALL ctl_stop( 'STOP', 'Charnock based wind stress can be used in coupled mode only' ) 426 ENDIF 427 ALLOCATE( charn(jpi,jpj) ) 428 charn(:,:) = 0._wp 429 ENDIF 430 IF( ln_taw ) THEN ! wind stress 431 IF( .NOT. cpl_taw ) THEN 432 CALL ctl_stop( 'STOP', 'wind stress from wave model can be used in coupled mode only, use ln_cdgw instead' ) 433 ENDIF 434 ALLOCATE( tawx(jpi,jpj) ) 435 ALLOCATE( tawy(jpi,jpj) ) 436 ALLOCATE( twox(jpi,jpj) ) 437 ALLOCATE( twoy(jpi,jpj) ) 438 ALLOCATE( tauoc_wavex(jpi,jpj) ) 439 ALLOCATE( tauoc_wavey(jpi,jpj) ) 440 tawx(:,:) = 0._wp 441 tawy(:,:) = 0._wp 442 twox(:,:) = 0._wp 443 twoy(:,:) = 0._wp 444 tauoc_wavex(:,:) = 1._wp 445 tauoc_wavey(:,:) = 1._wp 446 ENDIF 447 448 IF( ln_phioc ) THEN ! TKE flux 449 IF( .NOT. cpl_phioc ) THEN 450 CALL ctl_stop( 'STOP', 'phioc can be used in coupled mode only' ) 451 ENDIF 452 ALLOCATE( phioc(jpi,jpj) ) 453 phioc(:,:) = 0._wp 454 ENDIF 455 456 IF( ln_tauoc ) THEN ! normalized wave stress into the ocean 457 IF( .NOT. cpl_wstrf ) THEN 458 ALLOCATE( sf_tauoc(1), STAT=ierror ) !* allocate and fill sf_wave with sn_tauoc 459 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_tauoc structure' ) 460 ! 461 ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1) ) 462 IF( sn_tauoc%ln_tint ) ALLOCATE( sf_tauoc(1)%fdta(jpi,jpj,1,2) ) 463 CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave_init', 'Wave module', 'namsbc_wave' ) 464 ENDIF 465 ALLOCATE( tauoc_wave(jpi,jpj) ) 466 tauoc_wave(:,:) = 0._wp 467 ENDIF 468 469 IF( ln_sdw ) THEN ! Stokes drift 470 ! 1. Find out how many fields have to be read from file if not coupled 471 jpfld=0 472 jp_usd=0 ; jp_vsd=0 ; jp_hsw=0 ; jp_wmp=0 473 IF( .NOT. cpl_sdrftx ) THEN 474 jpfld = jpfld + 1 475 jp_usd = jpfld 476 ENDIF 477 IF( .NOT. cpl_sdrfty ) THEN 478 jpfld = jpfld + 1 479 jp_vsd = jpfld 480 ENDIF 481 IF( .NOT. cpl_hsig ) THEN 482 jpfld = jpfld + 1 483 jp_hsw = jpfld 484 ENDIF 485 IF( .NOT. cpl_wper ) THEN 486 jpfld = jpfld + 1 487 jp_wmp = jpfld 488 ENDIF 489 ! 2. Read from file only the non-coupled fields 490 IF( jpfld > 0 ) THEN 491 ALLOCATE( slf_i(jpfld) ) 492 IF( jp_usd > 0 ) slf_i(jp_usd) = sn_usd 493 IF( jp_vsd > 0 ) slf_i(jp_vsd) = sn_vsd 494 IF( jp_hsw > 0 ) slf_i(jp_hsw) = sn_hsw 495 IF( jp_wmp > 0 ) slf_i(jp_wmp) = sn_wmp 496 ALLOCATE( sf_sd(jpfld), STAT=ierror ) !* allocate and fill sf_sd with stokes drift 497 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' ) 498 ! 499 DO ifpr= 1, jpfld 500 ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) ) 501 IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) ) 502 END DO 503 ! 504 CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave_init', 'Wave module ', 'namsbc_wave' ) 505 ENDIF 386 506 ! 387 ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1) ) 388 IF( sn_cdg%ln_tint ) ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) ) 389 CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave_init', 'Wave module ', 'namsbc_wave' ) 390 ENDIF 391 ALLOCATE( cdn_wave(jpi,jpj) ) 392 ENDIF 393 394 IF( ln_tauwoc ) THEN 395 IF( .NOT. cpl_tauwoc ) THEN 396 ALLOCATE( sf_tauwoc(1), STAT=ierror ) !* allocate and fill sf_wave with sn_tauwoc 397 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' ) 507 ! 3. Wave number (only needed for Qiao parametrisation, ln_zdfqiao=T) 508 IF( .NOT. cpl_wnum ) THEN 509 ALLOCATE( sf_wn(1), STAT=ierror ) !* allocate and fill sf_wave with sn_wnum 510 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wn structure' ) 511 ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1) ) 512 IF( sn_wnum%ln_tint ) ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) ) 513 CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' ) 514 ENDIF 398 515 ! 399 ALLOCATE( sf_tauwoc(1)%fnow(jpi,jpj,1) ) 400 IF( sn_tauwoc%ln_tint ) ALLOCATE( sf_tauwoc(1)%fdta(jpi,jpj,1,2) ) 401 CALL fld_fill( sf_tauwoc, (/ sn_tauwoc /), cn_dir, 'sbc_wave_init', 'Wave module', 'namsbc_wave' ) 402 ENDIF 403 ALLOCATE( tauoc_wave(jpi,jpj) ) 404 ENDIF 405 406 IF( ln_tauw ) THEN 407 IF( .NOT. cpl_tauw ) THEN 408 ALLOCATE( sf_tauw(2), STAT=ierror ) !* allocate and fill sf_wave with sn_tauwx/y 409 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_tauw structure' ) 410 ! 411 ALLOCATE( slf_j(2) ) 412 slf_j(1) = sn_tauwx 413 slf_j(2) = sn_tauwy 414 ALLOCATE( sf_tauw(1)%fnow(jpi,jpj,1) ) 415 ALLOCATE( sf_tauw(2)%fnow(jpi,jpj,1) ) 416 IF( slf_j(1)%ln_tint ) ALLOCATE( sf_tauw(1)%fdta(jpi,jpj,1,2) ) 417 IF( slf_j(2)%ln_tint ) ALLOCATE( sf_tauw(2)%fdta(jpi,jpj,1,2) ) 418 CALL fld_fill( sf_tauw, (/ slf_j /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' ) 419 ENDIF 420 ALLOCATE( tauw_x(jpi,jpj) ) 421 ALLOCATE( tauw_y(jpi,jpj) ) 422 ENDIF 423 424 IF( ln_sdw ) THEN ! Find out how many fields have to be read from file if not coupled 425 jpfld=0 426 jp_usd=0 ; jp_vsd=0 ; jp_hsw=0 ; jp_wmp=0 ; jp_wfr=0 427 IF( .NOT. cpl_sdrftx ) THEN 428 jpfld = jpfld + 1 429 jp_usd = jpfld 430 ENDIF 431 IF( .NOT. cpl_sdrfty ) THEN 432 jpfld = jpfld + 1 433 jp_vsd = jpfld 434 ENDIF 435 IF( .NOT. cpl_hsig .AND. ll_st_bv_li ) THEN 436 jpfld = jpfld + 1 437 jp_hsw = jpfld 438 ENDIF 439 IF( .NOT. cpl_wper .AND. ll_st_bv_li ) THEN 440 jpfld = jpfld + 1 441 jp_wmp = jpfld 442 ENDIF 443 IF( .NOT. cpl_wfreq .AND. ll_st_peakfr ) THEN 444 jpfld = jpfld + 1 445 jp_wfr = jpfld 446 ENDIF 447 448 ! Read from file only the non-coupled fields 449 IF( jpfld > 0 ) THEN 450 ALLOCATE( slf_i(jpfld) ) 451 IF( jp_usd > 0 ) slf_i(jp_usd) = sn_usd 452 IF( jp_vsd > 0 ) slf_i(jp_vsd) = sn_vsd 453 IF( jp_hsw > 0 ) slf_i(jp_hsw) = sn_hsw 454 IF( jp_wmp > 0 ) slf_i(jp_wmp) = sn_wmp 455 IF( jp_wfr > 0 ) slf_i(jp_wfr) = sn_wfr 456 457 ALLOCATE( sf_sd(jpfld), STAT=ierror ) !* allocate and fill sf_sd with stokes drift 458 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' ) 459 ! 460 DO ifpr= 1, jpfld 461 ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) ) 462 IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) ) 463 END DO 464 ! 465 CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave_init', 'Wave module ', 'namsbc_wave' ) 466 ENDIF 467 ALLOCATE( usd (jpi,jpj,jpk), vsd (jpi,jpj,jpk), wsd(jpi,jpj,jpk) ) 468 ALLOCATE( hsw (jpi,jpj) , wmp (jpi,jpj) ) 469 ALLOCATE( wfreq(jpi,jpj) ) 470 ALLOCATE( ut0sd(jpi,jpj) , vt0sd(jpi,jpj) ) 471 ALLOCATE( div_sd(jpi,jpj) ) 472 ALLOCATE( tsd2d (jpi,jpj) ) 473 474 ut0sd(:,:) = 0._wp 475 vt0sd(:,:) = 0._wp 476 hsw(:,:) = 0._wp 477 wmp(:,:) = 0._wp 478 479 usd(:,:,:) = 0._wp 480 vsd(:,:,:) = 0._wp 481 wsd(:,:,:) = 0._wp 482 ! Wave number needed only if ln_zdfswm=T 483 IF( .NOT. cpl_wnum ) THEN 484 ALLOCATE( sf_wn(1), STAT=ierror ) !* allocate and fill sf_wave with sn_wnum 485 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave_init: unable toallocate sf_wave structure' ) 486 ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1) ) 487 IF( sn_wnum%ln_tint ) ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) ) 488 CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' ) 489 ENDIF 490 ALLOCATE( wnum(jpi,jpj) ) 516 ENDIF 517 ! 491 518 ENDIF 492 519 ! -
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/ZDF/zdfsh2.F90
r12377 r12991 13 13 USE oce 14 14 USE dom_oce ! domain: ocean 15 USE sbcwave ! Surface Waves (add Stokes shear) 16 USE sbc_oce , ONLY: ln_stshear !Stoked Drift shear contribution 15 17 ! 16 18 USE in_out_manager ! I/O manager … … 21 23 22 24 PUBLIC zdf_sh2 ! called by zdftke, zdfglf, and zdfric 23 25 24 26 !! * Substitutions 25 27 # include "do_loop_substitute.h90" … … 58 60 !!-------------------------------------------------------------------- 59 61 ! 60 DO jk = 2, jpkm1 61 DO_2D_10_10 62 zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) ) & 63 & * ( uu(ji,jj,jk-1,Kmm) - uu(ji,jj,jk,Kmm) ) & 64 & * ( uu(ji,jj,jk-1,Kbb) - uu(ji,jj,jk,Kbb) ) / ( e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) ) * wumask(ji,jj,jk) 65 zsh2v(ji,jj) = ( p_avm(ji,jj+1,jk) + p_avm(ji,jj,jk) ) & 66 & * ( vv(ji,jj,jk-1,Kmm) - vv(ji,jj,jk,Kmm) ) & 67 & * ( vv(ji,jj,jk-1,Kbb) - vv(ji,jj,jk,Kbb) ) / ( e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) ) * wvmask(ji,jj,jk) 68 END_2D 62 DO jk = 2, jpkm1 !* Shear production at uw- and vw-points (energy conserving form) 63 IF ( cpl_sdrftx .AND. ln_stshear ) THEN ! Surface Stokes Drift available ===>>> shear + stokes drift contibution 64 DO_2D_10_10 65 zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) ) & 66 & * ( uu (ji,jj,jk-1,Kmm) - uu (ji,jj,jk,Kmm) & 67 & + usd(ji,jj,jk-1) - usd(ji,jj,jk) ) & 68 & * ( uu (ji,jj,jk-1,Kbb) - uu (ji,jj,jk,Kbb) ) & 69 & / ( e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) ) * wumask(ji,jj,jk) 70 zsh2v(ji,jj) = ( p_avm(ji,jj+1,jk) + p_avm(ji,jj,jk) ) & 71 & * ( vv (ji,jj,jk-1,Kmm) - vv (ji,jj,jk,Kmm) & 72 & + vsd(ji,jj,jk-1) - vsd(ji,jj,jk) ) & 73 & * ( vv (ji,jj,jk-1,Kbb) - vv (ji,jj,jk,Kbb) ) & 74 &/ ( e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) ) * wvmask(ji,jj,jk) 75 END_2D 76 ELSE 77 DO_2D_10_10 78 zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) ) & 79 & * ( uu(ji,jj,jk-1,Kmm) - uu(ji,jj,jk,Kmm) ) & 80 & * ( uu(ji,jj,jk-1,Kbb) - uu(ji,jj,jk,Kbb) ) & 81 & / ( e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) ) * wumask(ji,jj,jk) 82 zsh2v(ji,jj) = ( p_avm(ji,jj+1,jk) + p_avm(ji,jj,jk) ) & 83 & * ( vv(ji,jj,jk-1,Kmm) - vv(ji,jj,jk,Kmm) ) & 84 & * ( vv(ji,jj,jk-1,Kbb) - vv(ji,jj,jk,Kbb) ) & 85 & / ( e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) ) * wvmask(ji,jj,jk) 86 END_2D 87 ENDIF 69 88 DO_2D_00_00 70 89 p_sh2(ji,jj,jk) = 0.25 * ( ( zsh2u(ji-1,jj) + zsh2u(ji,jj) ) * ( 2. - umask(ji-1,jj,jk) * umask(ji,jj,jk) ) & 71 90 & + ( zsh2v(ji,jj-1) + zsh2v(ji,jj) ) * ( 2. - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) ) ) 72 91 END_2D 73 END DO 92 END DO 74 93 ! 75 94 END SUBROUTINE zdf_sh2 -
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/ZDF/zdftke.F90
r12702 r12991 52 52 USE prtctl ! Print control 53 53 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 54 USE sbcwave ! Surface boundary waves 54 55 55 56 IMPLICIT NONE … … 62 63 ! !!** Namelist namzdf_tke ** 63 64 LOGICAL :: ln_mxl0 ! mixing length scale surface value as function of wind stress or not 65 LOGICAL :: ln_mxhsw ! mixing length scale surface value as a fonction of wave height 64 66 INTEGER :: nn_mxl ! type of mixing length (=0/1/2/3) 65 67 REAL(wp) :: rn_mxl0 ! surface min value of mixing length (kappa*z_o=0.4*0.1 m) [m] … … 74 76 INTEGER :: nn_etau ! type of depth penetration of surface tke (=0/1/2/3) 75 77 INTEGER :: nn_htau ! type of tke profile of penetration (=0/1) 78 INTEGER :: nn_bc_surf! surface condition (0/1=Dir/Neum) ! Only applicable for wave coupling 79 INTEGER :: nn_bc_bot ! surface condition (0/1=Dir/Neum) ! Only applicable for wave coupling 76 80 REAL(wp) :: rn_efr ! fraction of TKE surface value which penetrates in the ocean 77 81 REAL(wp) :: rn_eice ! =0 ON below sea-ice, =4 OFF when ice fraction > 1/4 … … 201 205 REAL(wp) :: zus , zwlc , zind ! - - 202 206 REAL(wp) :: zzd_up, zzd_lw ! - - 207 REAL(wp) :: ztaui, ztauj, z1_norm 203 208 INTEGER , DIMENSION(jpi,jpj) :: imlc 204 REAL(wp), DIMENSION(jpi,jpj) :: zhlc, zfr_i 209 REAL(wp), DIMENSION(jpi,jpj) :: zhlc, zfr_i, zWlc2 205 210 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpelc, zdiag, zd_up, zd_lw 206 211 !!-------------------------------------------------------------------- … … 210 215 zfact2 = 1.5_wp * rn_Dt * rn_ediss 211 216 zfact3 = 0.5_wp * rn_ediss 217 ! 218 zpelc(:,:,:) = 0._wp ! need to be initialised in case ln_lc is not used 212 219 ! 213 220 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 253 260 ! 254 261 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 255 IF( ln_lc ) THEN ! Langmuir circulation source term added to tke !(Axell JGR 2002)262 IF( ln_lc ) THEN ! Langmuir circulation source term added to tke (Axell JGR 2002) 256 263 ! !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 257 264 ! 258 ! !* total energy produce by LC : cumulative sum over jk 265 ! !* Langmuir velocity scale 266 ! 267 IF ( cpl_sdrftx ) THEN ! Surface Stokes Drift available 268 ! ! Craik-Leibovich velocity scale Wlc = ( u* u_s )^1/2 with u* = (taum/rho0)^1/2 269 ! ! associated kinetic energy : 1/2 (Wlc)^2 = u* u_s 270 ! ! more precisely, it is the dot product that must be used : 271 ! ! 1/2 (W_lc)^2 = MAX( u* u_s + v* v_s , 0 ) only the positive part 272 !!gm ! PS: currently we don't have neither the 2 stress components at t-point !nor the angle between u* and u_s 273 !!gm ! so we will overestimate the LC velocity.... !!gm I will do the work if !LC have an effect ! 274 DO_2D_00_00 275 !!XC zWlc2(ji,jj) = 0.5_wp * SQRT( taum(ji,jj) * r1_rho0 * ( ut0sd(ji,jj)**2 +vt0sd(ji,jj)**2 ) ) 276 zWlc2(ji,jj) = 0.5_wp * ( ut0sd(ji,jj)**2 +vt0sd(ji,jj)**2 ) 277 END_2D 278 ! 279 ! Projection of Stokes drift in the wind stress direction 280 ! 281 DO_2D_00_00 282 ztaui = 0.5_wp * ( utau(ji,jj) + utau(ji-1,jj) ) 283 ztauj = 0.5_wp * ( vtau(ji,jj) + vtau(ji,jj-1) ) 284 z1_norm = 1._wp / MAX( SQRT(ztaui*ztaui+ztauj*ztauj), 1.e-12 ) * tmask(ji,jj,1) 285 zWlc2(ji,jj) = 0.5_wp * z1_norm * ( MAX( ut0sd(ji,jj)*ztaui + vt0sd(ji,jj)*ztauj, 0._wp ) )**2 286 END_2D 287 CALL lbc_lnk ( 'zdftke', zWlc2, 'T', 1. ) 288 ! 289 ELSE ! Surface Stokes drift deduced from surface stress 290 ! ! Wlc = u_s with u_s = 0.016*U_10m, the surface stokes drift (Axell 2002, Eq.44) 291 ! ! using |tau| = rho_air Cd |U_10m|^2 , it comes: 292 ! ! Wlc = 0.016 * [|tau|/(rho_air Cdrag) ]^1/2 and thus: 293 ! ! 1/2 Wlc^2 = 0.5 * 0.016 * 0.016 |tau| /( rho_air Cdrag ) 294 zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag ) ! to convert stress in 10m wind using a constant drag 295 DO_2D_11_11 296 zWlc2(ji,jj) = zcof * taum(ji,jj) 297 END_2D 298 ! 299 ENDIF 300 ! 301 ! !* Depth of the LC circulation (Axell 2002, Eq.47) 302 ! !- LHS of Eq.47 259 303 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0._wp ) * gdepw(:,:,1,Kmm) * e3w(:,:,1,Kmm) 260 304 DO jk = 2, jpk 261 305 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * gdepw(:,:,jk,Kmm) * e3w(:,:,jk,Kmm) 262 306 END DO 263 ! !* finite Langmuir Circulation depth264 zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )307 ! 308 ! !- compare LHS to RHS of Eq.47 265 309 imlc(:,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point (=2 over land) 266 310 DO_3DS_11_11( jpkm1, 2, -1 ) 267 zus = zcof * taum(ji,jj) 268 IF( zpelc(ji,jj,jk) > zus ) imlc(ji,jj) = jk 311 IF( zpelc(ji,jj,jk) > zWlc2(ji,jj) ) imlc(ji,jj) = jk 269 312 END_3D 270 313 ! ! finite LC depth … … 272 315 zhlc(ji,jj) = gdepw(ji,jj,imlc(ji,jj),Kmm) 273 316 END_2D 317 ! 274 318 zcof = 0.016 / SQRT( zrhoa * zcdrag ) 275 319 DO_2D_00_00 276 zus = zcof * SQRT( taum(ji,jj) )! Stokes drift320 zus = SQRT( 2. * zWlc2(ji,jj) ) ! Stokes drift 277 321 zfr_i(ji,jj) = ( 1._wp - 4._wp * fr_i(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 278 322 IF (zfr_i(ji,jj) < 0. ) zfr_i(ji,jj) = 0. … … 327 371 & ) * wmask(ji,jj,jk) 328 372 END_3D 373 ! 374 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 375 ! ! choose to keep coherence with previous estimation of 376 ! ! Surface boundary condition on tke 377 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 378 ! [EC] Should we keep this?? 379 IF ( ln_isfcav ) THEN 380 DO_2D_11_11 ! en(mikt(ji,jj)) = rn_emin 381 en(ji,jj,mikt(ji,jj))=rn_emin * tmask(ji,jj,1) 382 END_2D 383 END IF 384 385 IF ( cpl_phioc .and. ln_phioc ) THEN 386 SELECT CASE (nn_bc_surf) !! Dirichlet Boundary Condition using surface TKE flux from waves 387 388 CASE ( 0 ) 389 390 DO_2D_00_00 ! en(1) = rn_ebb taum / rho0 (min value rn_emin0) 391 IF ( phioc(ji,jj) < 0 ) phioc(ji,jj) = 0._wp 392 en(ji,jj,1) = MAX( rn_emin0, .5 * ( 15.8 * phioc(ji,jj) / rho0 )**(2./3.) ) * tmask(ji,jj,1) 393 zdiag(ji,jj,1) = 1._wp/en(ji,jj,1) ! choose to keep coherence with former estimation of 394 zd_lw(ji,jj,1) = 1._wp ! zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) 395 zd_up(ji,jj,1) = 0._wp 396 END_2D 397 398 CASE ( 1 ) 399 DO_2D_00_00 400 IF ( phioc(ji,jj) < 0 ) phioc(ji,jj) = 0._wp 401 en(ji,jj,2) = en(ji,jj,2) + ( rn_Dt * phioc(ji,jj) / rho0 ) /e3w(ji,jj,2,Kmm) 402 en(ji,jj,1) = en(ji,jj,2) + (2 * e3t(ji,jj,1,Kmm) * phioc(ji,jj)) / ( p_avm(ji,jj,1) + p_avm(ji,jj,2) ) 403 zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2) 404 zd_lw(ji,jj,2) = 0._wp 405 zd_up(ji,jj,1) = 0._wp 406 END_2D 407 ! 408 END SELECT 409 410 ELSE ! TKE Dirichlet boundary condition (without wave coupling) 411 412 DO_2D_00_00 ! en(1) = rn_ebb taum / rho0 (min value rn_emin0) 413 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 414 zdiag(ji,jj,1) = 1._wp/en(ji,jj,1) ! choose to keep coherence with former estimation of 415 zd_lw(ji,jj,1) = 1._wp ! zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) 416 zd_up(ji,jj,1) = 0._wp 417 END_2D 418 419 ENDIF 420 ! 329 421 ! !* Matrix inversion from level 2 (tke prescribed at level 1) 330 DO_3D_00_00( 3, jpkm1 ) 422 ! DO_3D_00_00( 3, jpkm1 ) 423 DO_3D_00_00( 2, jpkm1 ) ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 331 424 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 332 425 END_3D 333 DO_2D_00_00 334 zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke 335 END_2D 336 DO_3D_00_00( 3, jpkm1 ) 426 !XC : commented to allow for neumann boundary condition 427 ! DO_2D_00_00 428 ! zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke 429 ! END_2D 430 ! DO_3D_00_00( 3, jpkm1 ) 431 DO_3D_00_00( 2, jpkm1 ) ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 337 432 zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1) 338 433 END_3D … … 340 435 en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 341 436 END_2D 342 DO_3DS_00_00( jpk-2, 2, -1 ) 437 DO_3DS_00_00( jpk-2, 2, -1 ) ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 343 438 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 344 439 END_3D 345 DO_3D_00_00( 2, jpkm1 ) 440 DO_3D_00_00( 2, jpkm1 ) ! set the minimum value of tke 346 441 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) 347 442 END_3D … … 436 531 zmxld(:,:,:) = rmxl_min 437 532 ! 438 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g) 439 zraug = vkarmn * 2.e5_wp / ( rho0 * grav ) 440 DO_2D_00_00 441 zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) ) 442 END_2D 443 ELSE 444 zmxlm(:,:,1) = rn_mxl0 533 IF(ln_sdw .AND. ln_mxhsw) THEN 534 zmxlm(:,:,1)= vkarmn * MAX ( 1.6 * hsw(:,:) , 0.02 ) ! surface mixing length = F(wave height) 535 ! from terray et al 1999 and mellor and blumberg 2004 it should be 0.85 and not 1.6 536 zcoef = vkarmn * ( (rn_ediff*rn_ediss)**0.25 ) / rn_ediff 537 zmxlm(:,:,1)= zcoef * MAX ( 1.6 * hsw(:,:) , 0.02 ) ! surface mixing length = F(wave height) 538 ELSE 539 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g) 540 zraug = vkarmn * 2.e5_wp / ( rho0 * grav ) 541 DO_2D_00_00 542 zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) ) 543 END_2D 544 ELSE 545 zmxlm(:,:,1) = rn_mxl0 546 ENDIF 445 547 ENDIF 446 548 ! … … 550 652 & rn_emin0, rn_bshear, nn_mxl , ln_mxl0 , & 551 653 & rn_mxl0 , nn_pdl , ln_drg , ln_lc , rn_lc, & 552 & nn_etau , nn_htau , rn_efr , rn_eice 654 & nn_etau , nn_htau , rn_efr , rn_eice , & 655 & nn_bc_surf, nn_bc_bot, ln_mxhsw 553 656 !!---------------------------------------------------------------------- 554 657 !
Note: See TracChangeset
for help on using the changeset viewer.