Changeset 7807
- Timestamp:
- 2017-03-17T10:44:05+01:00 (7 years ago)
- Location:
- branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 1 added
- 11 edited
-
BDY/bdyini.F90 (modified) (1 diff)
-
DIA/diawri.F90 (modified) (1 diff)
-
DYN/dynstcor.F90 (added)
-
SBC/sbc_oce.F90 (modified) (3 diffs)
-
SBC/sbcblk_core.F90 (modified) (10 diffs)
-
SBC/sbcblk_mfs.F90 (modified) (1 diff)
-
SBC/sbcflx.F90 (modified) (3 diffs)
-
SBC/sbcmod.F90 (modified) (6 diffs)
-
SBC/sbcwave.F90 (modified) (4 diffs)
-
ZDF/zdfgls.F90 (modified) (14 diffs)
-
step.F90 (modified) (1 diff)
-
step_oce.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90
r6204 r7807 297 297 IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.' 298 298 ln_tra_dmp(ib_bdy)=.false. 299 ELSEIF ( cn_tra(ib_bdy) == 'frs' ) THEN300 CALL ctl_stop( 'Use FRS OR relaxation' )301 299 ELSE 302 300 IF(lwp) WRITE(numout,*) ' + T/S relaxation zone' -
branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
r6204 r7807 767 767 768 768 IF( lk_vvl ) THEN 769 CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * fse3t_n(:,:,:), ndim_T , ndex_T ) ! heat content770 CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * fse3t_n(:,:,:), ndim_T , ndex_T ) ! salt content769 CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem), ndim_T , ndex_T ) ! heat content 770 CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal), ndim_T , ndex_T ) ! salt content 771 771 CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * fse3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface heat content 772 772 CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * fse3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface salinity content -
branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
r5407 r7807 62 62 ! !: = 1 global mean of e-p-r set to zero at each nn_fsbc time step 63 63 ! !: = 2 annual global mean of e-p-r set to zero 64 LOGICAL , PUBLIC :: ln_wave !: true if some coupling with wave model 65 LOGICAL , PUBLIC :: ln_cdgw !: true if neutral drag coefficient from wave model 66 LOGICAL , PUBLIC :: ln_sdw !: true if 3d stokes drift from wave model 64 ! LOGICAL , PUBLIC :: ln_wave !: true if some coupling with wave model 65 LOGICAL , PUBLIC :: ln_cdgw = .FALSE. !: true if neutral drag coefficient from wave model 66 ! LOGICAL , PUBLIC :: ln_sdw !: true if 3d stokes drift from wave model 67 ! WAVE2NEMO 11.12.2016 68 LOGICAL , PUBLIC :: ln_wavetke = .FALSE. !: true if wave parametersare read 69 LOGICAL , PUBLIC :: ln_stcor = .FALSE. !: true if Stokes-Coriolisforcing is included 70 LOGICAL , PUBLIC :: ln_tauoc = .FALSE. !: true if wave-modifiedwater-side stress is used 71 ! END WAVE2NEMO 72 73 67 74 ! 68 75 LOGICAL , PUBLIC :: ln_icebergs !: Icebergs … … 125 132 #endif 126 133 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xcplmask !: coupling mask for ln_mixcpl (warning: allocated in sbccpl) 127 134 ! START WAVE2NEMO 26.11.2016 135 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tauoc_wavepar !:normalized stress into the ocean from wave model CCC 136 REAL(wp), PUBLIC,ALLOCATABLE,SAVE,DIMENSION (:,:) :: wspd_wavepar !CCCmoved to sbc_oce 137 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: cdn_wave !:wave model drag coefficient [N/m2] 138 ! END 139 !WAVE2NEMO (14.12.2016) rn_crban moved to here, it needs to be a more global variable 140 REAL(wp),PUBLIC,ALLOCATABLE,SAVE, DIMENSION (:,:) :: rn_crban 128 141 !!---------------------------------------------------------------------- 129 142 !! Sea Surface Mean fields … … 176 189 ALLOCATE( e3t_m(jpi,jpj) , STAT=ierr(5) ) 177 190 #endif 178 ! 191 ! Start WAVE2NEMO 11.12.2016 192 ALLOCATE( cdn_wave(jpi,jpj) ) 193 ALLOCATE( wspd_wavepar(jpi,jpj) ) 194 ! End 195 ! 196 ALLOCATE( rn_crban(jpi,jpj) ) !WAVE2NEMO (13.12.2016) 197 198 ! 179 199 sbc_oce_alloc = MAXVAL( ierr ) 180 200 IF( lk_mpp ) CALL mpp_sum ( sbc_oce_alloc ) -
branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
r5582 r7807 17 17 !! 3.7 ! 2014-06 (L. Brodeau) simplification and optimization of CORE bulk 18 18 !!---------------------------------------------------------------------- 19 19 !! 3.6 2016-06 (V.Alari) Added wave dependent drag 20 !coefficient and ocean side stress 20 21 !!---------------------------------------------------------------------- 21 22 !! sbc_blk_core : bulk formulation as ocean surface boundary condition (forced mode, CORE bulk formulea) … … 41 42 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 42 43 USE prtctl ! Print control 43 USE sbcwave, ONLY : cdn_wave ! wave module 44 ! USE sbcwave, ONLY : cdn_wave ! wave module (V.Alari 27.06.2016) We get 45 ! info from sbc_oce 44 46 USE sbc_ice ! Surface boundary condition: ice fields 45 47 USE lib_fortran ! to use key_nosignedzero … … 183 185 ! 184 186 lhftau = ln_taudif ! do we use HF tau information? 187 188 189 185 190 jfld = jpfld - COUNT( (/.NOT. lhftau/) ) 186 191 ! … … 192 197 END DO 193 198 ! ! fill sf with slf_i and control print 199 194 200 CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_core', 'flux formulation for ocean surface boundary condition', 'namsbc_core' ) 195 201 ! … … 258 264 REAL(wp), DIMENSION(:,:), POINTER :: zq_zu ! air spec. hum. at wind speed height 259 265 !!--------------------------------------------------------------------- 266 REAL(wp) :: ztheta ! local variable, wind direction,V.Alari 27.06.2016 267 260 268 ! 261 269 IF( nn_timing == 1 ) CALL timing_start('blk_oce_core') … … 285 293 END DO 286 294 #endif 295 ! Start 27.06.2016 (V.Alari) 296 ! Wind vector at T points 297 IF ( ln_cdgw ) THEN 298 ! Use neutral 10-m wind speed from wave model if available 299 DO jj = 2, jpjm1 300 DO ji = fs_2, fs_jpim1 ! vect. opt. 301 ! Local wind direction from non-neutral 10-m wind vector 302 ! (relative to grid and no current) 303 ! since we do not have the directional information from the wave 304 ! model 305 ztheta = ATAN2(sf(jp_wndi)%fnow(ji,jj,1),sf(jp_wndj)%fnow(ji,jj,1)) 306 ! Wind vector magnitude from 10-m neutral wind speed from wave 307 ! model 308 zwnd_i(ji,jj) = wspd_wavepar(ji,jj) * SIN(ztheta) 309 zwnd_j(ji,jj) = wspd_wavepar(ji,jj) * COS(ztheta) 310 ! Correct for surface current, 0.0 <= rn_rrelwind <= 1.0 311 zwnd_i(ji,jj) = (zwnd_i(ji,jj) - 0.5*rn_vfac*(pu(ji-1,jj )+pu(ji,jj))) 312 zwnd_j(ji,jj) = (zwnd_j(ji,jj) - 0.5*rn_vfac*(pv(ji ,jj-1)+pv(ji,jj))) 313 END DO 314 END DO 315 ELSE ! If no neutral wind from wave model 287 316 DO jj = 2, jpjm1 288 317 DO ji = fs_2, fs_jpim1 ! vect. opt. … … 291 320 END DO 292 321 END DO 322 ENDIF ! End (V.Alari) 323 293 324 CALL lbc_lnk( zwnd_i(:,:) , 'T', -1. ) 294 325 CALL lbc_lnk( zwnd_j(:,:) , 'T', -1. ) … … 320 351 321 352 ! ... tau module, i and j component 353 ! Start 27.06.2016 (V.Alari) 354 IF ( ln_cdgw ) THEN 355 DO jj = 1, jpj 356 DO ji = 1, jpi 357 !zztmp = rhoa * wspd_wavepar(ji,jj) * Cd_n10(ji,jj) 358 zztmp = rhoa * wndm(ji,jj) * cdn_wave(ji,jj) 359 !taum(ji,jj) = tauoc_wavepar(ji,jj) 360 !taum(ji,jj) = zztmp * wspd_wavepar(ji,jj) 361 taum (ji,jj) = zztmp * wndm (ji,jj) 362 zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 363 zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 364 ENDDO 365 ENDDO 366 ELSE 322 367 DO jj = 1, jpj 323 368 DO ji = 1, jpi … … 328 373 END DO 329 374 END DO 330 375 ENDIF ! End, (V.Alari) 331 376 ! ... add the HF tau contribution to the wind stress module? 332 377 IF( lhftau ) THEN … … 349 394 CALL lbc_lnk( vtau(:,:), 'V', -1. ) 350 395 351 396 ! START: 27.06.2015 (V. Alari). Added tauoc 397 398 IF (ln_tauoc) THEN 399 utau(:,:) = utau(:,:)*tauoc_wavepar(:,:) 400 vtau(:,:) = vtau(:,:)*tauoc_wavepar(:,:) 401 taum(:,:) = taum(:,:)*tauoc_wavepar(:,:) 402 ENDIF 403 404 ! END (V.Alari) 352 405 ! Turbulent fluxes over ocean 353 406 ! ----------------------------- -
branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90
r5215 r7807 24 24 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 25 25 USE prtctl ! Print control 26 USE sbcwave,ONLY : cdn_wave !wave module26 ! USE sbcwave,ONLY : cdn_wave !wave module 27 27 28 28 IMPLICIT NONE -
branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90
r4990 r7807 22 22 USE lib_mpp ! distribued memory computing library 23 23 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 24 USE sbcwave ! Wave parameters 24 25 25 26 IMPLICIT NONE … … 152 153 zty = vtau(ji ,jj-1) + vtau(ji,jj) 153 154 zmod = 0.5 * SQRT( ztx * ztx + zty * zty ) 155 IF (ln_cdgw) THEN 156 zcdrag = 1.5e-3 157 IF (cdn_wave(ji,jj)>=0.) zcdrag = cdn_wave(ji,jj) 158 zcoef = 1./(zrhoa*zcdrag) 159 ENDIF 154 160 taum(ji,jj) = zmod 155 161 wndm(ji,jj) = SQRT( zmod * zcoef ) … … 158 164 taum(:,:) = taum(:,:) * tmask(:,:,1) ; wndm(:,:) = wndm(:,:) * tmask(:,:,1) 159 165 CALL lbc_lnk( taum(:,:), 'T', 1. ) ; CALL lbc_lnk( wndm(:,:), 'T', 1. ) 166 167 IF (ln_tauoc) THEN 168 utau(:,:) = utau(:,:)*tauoc_wavepar(:,:) 169 vtau(:,:) = vtau(:,:)*tauoc_wavepar(:,:) 170 taum(:,:) = taum(:,:)*tauoc_wavepar(:,:) 171 ENDIF 160 172 161 173 IF( nitend-nit000 <= 100 .AND. lwp ) THEN ! control print (if less than 100 time-step asked) -
branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
r5628 r7807 85 85 !!---------------------------------------------------------------------- 86 86 INTEGER :: icpt ! local integer 87 !! 87 !! V.Alari 26.06.2016 Added logicals for waves, deleted ln_wave ln_sdw, which came 88 !with UKMO 88 89 NAMELIST/namsbc/ nn_fsbc , ln_ana , ln_flx, ln_blk_clio, ln_blk_core, ln_mixcpl, & 89 90 & ln_blk_mfs, ln_apr_dyn, nn_ice, nn_ice_embd, ln_dm2dc , ln_rnf , & 90 & ln_ssr , nn_isf , nn_fwb, ln_cdgw , ln_wave , ln_sdw , & 91 & nn_lsm , nn_limflx , nn_components, ln_cpl 91 & ln_ssr , nn_isf , nn_fwb, & 92 & nn_lsm , nn_limflx , nn_components, ln_cpl,ln_stcor, ln_wavetke, & 93 & ln_tauoc, ln_cdgw 92 94 INTEGER :: ios 93 95 INTEGER :: ierr, ierr0, ierr1, ierr2, ierr3, jpm … … 120 122 nn_ice = 0 121 123 ENDIF 122 124 ! V.Alari added wave logicalc (26.06.2016) 123 125 IF(lwp) THEN ! Control print 124 126 WRITE(numout,*) ' Namelist namsbc (partly overwritten with CPP key setting)' … … 146 148 WRITE(numout,*) ' closed sea (=0/1) (set in namdom) nn_closea = ', nn_closea 147 149 WRITE(numout,*) ' n. of iterations if land-sea-mask applied nn_lsm = ', nn_lsm 148 ENDIF 150 WRITE(numout,*) ' Stokes-Coriolis forcing from wavemodel ln_stcor = ', ln_stcor 151 WRITE(numout,*) ' TKE wave breaking BC from wave model ln_wavetke = ', ln_wavetke 152 WRITE(numout,*) ' Wave-modified stress from wave model ln_tauoc = ', ln_tauoc 153 WRITE(numout,*) ' Drag coefficient from wave model ln_cdgw = ', ln_cdgw 154 155 156 ENDIF 149 157 150 158 ! LIM3 Multi-category heat flux formulation … … 179 187 180 188 ! ! Checks: 181 IF( nn_isf .EQ. 0 ) THEN ! variable initialisation if noice shelf189 IF( nn_isf .EQ. 0 ) THEN ! no specific treatment in vicinity of ice shelf 182 190 IF( sbc_isf_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_isf arrays' ) 183 191 fwfisf (:,:) = 0.0_wp ; fwfisf_b (:,:) = 0.0_wp … … 214 222 IF( ln_dm2dc .AND. .NOT.( ln_flx .OR. ln_blk_core ) .AND. nn_components /= jp_iam_opa ) & 215 223 & CALL ctl_stop( 'diurnal cycle into qsr field from daily values requires a flux or core-bulk formulation' ) 216 217 IF ( ln_wave ) THEN 224 225 ! V.Alari comment out UKMO wave stuff 226 ! IF ( ln_wave ) THEN 218 227 !Activated wave module but neither drag nor stokes drift activated 219 IF ( .NOT.(ln_cdgw .OR. ln_sdw) ) THEN220 CALL ctl_warn( 'Ask for wave coupling but nor drag coefficient (ln_cdgw=F) neither stokes drift activated (ln_sdw=F)' )228 ! IF ( .NOT.(ln_cdgw .OR. ln_sdw) ) THEN 229 ! CALL ctl_warn( 'Ask for wave coupling but nor drag coefficient (ln_cdgw=F) neither stokes drift activated (ln_sdw=F)' ) 221 230 !drag coefficient read from wave model definable only with mfs bulk formulae and core 222 ELSEIF (ln_cdgw .AND. .NOT.(ln_blk_mfs .OR. ln_blk_core) ) THEN 223 CALL ctl_stop( 'drag coefficient read from wave model definable only with mfs bulk formulae and core') 224 ENDIF 225 ELSE 226 IF ( ln_cdgw .OR. ln_sdw ) & 227 & CALL ctl_stop('Not Activated Wave Module (ln_wave=F) but & 228 & asked coupling with drag coefficient (ln_cdgw =T) or Stokes drift (ln_sdw=T) ') 229 ENDIF 231 ! ELSEIF (ln_cdgw .AND. .NOT.(ln_blk_mfs .OR. ln_blk_core) ) THEN 232 ! CALL ctl_stop( 'drag coefficient read from wave model definable only with mfs bulk formulae and core') 233 ! ENDIF 234 ! ELSE 235 ! IF ( ln_cdgw .OR. ln_sdw ) & 236 ! & CALL ctl_stop('Not Activated Wave Module (ln_wave=F) but & 237 ! & asked coupling with drag coefficient (ln_cdgw =T) or Stokes drift (ln_sdw=T) ') 238 ! ENDIF 239 ! End V.Alari 230 240 ! ! Choice of the Surface Boudary Condition (set nsbc) 231 241 ll_purecpl = ln_cpl .AND. .NOT. ln_mixcpl … … 352 362 IF( nn_components /= jp_iam_sas ) CALL sbc_ssm( kt ) ! ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m) 353 363 ! ! averaged over nf_sbc time-step 354 355 IF (ln_wave) CALL sbc_wave( kt ) 364 ! START V.Alari 26.06.2016 365 IF (ln_cdgw) CALL sbc_wave( kt ) 366 ! Read wave parameters if Stokes-Coriolis forcing OR wave parameters for 367 ! TKE are needed 368 ! Also read it if drag coefficient to ensure that we use the neutral 10m 369 ! from WAM. 370 IF (ln_stcor .OR. ln_wavetke .OR. ln_tauoc) THEN 371 CALL sbc_wavepar( kt ) 372 ENDIF 373 374 ! END 375 376 ! IF (ln_wave) CALL sbc_wave( kt ) 356 377 !== sbc formulation ==! 357 378 -
branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
r5215 r7807 4 4 !! Wave module 5 5 !!====================================================================== 6 !! History : 3.3.1 ! 2011-09 (Adani M) Original code: Drag Coefficient 7 !! : 3.4 ! 2012-10 (Adani M) Stokes Drift 6 !! History : 3.3.1 ! 2011-09 (Adani M) Original code 8 7 !!---------------------------------------------------------------------- 9 8 USE iom ! I/O manager library 10 9 USE in_out_manager ! I/O manager 11 10 USE lib_mpp ! distribued memory computing library 12 USE fldread ! read input fields 13 USE oce 14 USE sbc_oce ! Surface boundary condition: ocean fields 15 USE domvvl 16 17 11 USE fldread ! read input fields 12 USE sbc_oce ! Surface boundary condition: ocean fields 13 USE phycst ! physical constants 14 ! USE sbcblk_core, ONLY: Cd_n10 18 15 !!---------------------------------------------------------------------- 19 16 !! sbc_wave : read drag coefficient from wave model in netcdf files … … 23 20 PRIVATE 24 21 25 PUBLIC sbc_wave ! routine called in sbc_blk_core or sbc_blk_mfs 26 27 INTEGER , PARAMETER :: jpfld = 3 ! maximum number of files to read for srokes drift 28 INTEGER , PARAMETER :: jp_usd = 1 ! index of stokes drift (i-component) (m/s) at T-point 29 INTEGER , PARAMETER :: jp_vsd = 2 ! index of stokes drift (j-component) (m/s) at T-point 30 INTEGER , PARAMETER :: jp_wn = 3 ! index of wave number (1/m) at T-point 31 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_cd ! structure of input fields (file informations, fields read) Drag Coefficient 32 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_sd ! structure of input fields (file informations, fields read) Stokes Drift 33 REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:) :: cdn_wave 34 REAL(wp),ALLOCATABLE,DIMENSION (:,:) :: usd2d,vsd2d,uwavenum,vwavenum 35 REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:,:) :: usd3d,vsd3d,wsd3d 22 PUBLIC sbc_wavepar ! routine called in sbc_blk_core or sbc_blk_mfs 23 PUBLIC sbc_wave ! routine called in sbc_blk_core or sbc_blk_mfs 24 25 INTEGER , PARAMETER :: jpfld_wavepar = 7 ! maximum number of fields to be read 26 INTEGER , PARAMETER :: jp_ust = 1 ! index of Stokes velocity (east component) (m/s) at T-point 27 INTEGER , PARAMETER :: jp_vst = 2 ! index of Stokes velocity (north component) (m/s) at T-point 28 INTEGER , PARAMETER :: jp_swh = 3 ! index of significant wave height (m) 29 INTEGER , PARAMETER :: jp_mwp = 4 ! index of mean wave period (s) 30 INTEGER , PARAMETER :: jp_phioc = 5 ! index of normalized energy flux into the ocean (non-dim) 31 INTEGER , PARAMETER :: jp_tauoc = 6 ! index of normalized wave stress into the ocean (non-dim) 32 !INTEGER , PARAMETER :: jp_cdww = 7 ! index of wave drag coeff (non-dim) 33 INTEGER , PARAMETER :: jp_wspd = 7 ! index of 10 m neutral wind speed (m/s) 34 35 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_wave ! structure of input fields (file informations, fields read) 36 !REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:) :: cdn_wave ! CCC moved to sbc_oce 37 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_wavepar ! structure of input fields (file informations, fields read) 38 REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:) :: ust_wavepar, vst_wavepar, swh_wavepar, mwp_wavepar 39 REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:) :: phioc_wavepar !WAVE2NEMO 14.12.2016 moved rn_crban to sbc_oce 40 !REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:) :: tauoc_wavepar ! CCC moved to sbc_oce 41 !REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:) :: wspd_wavepar ! CCC moved to sbc_oce 42 !REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:) :: cdww_wavepar 36 43 37 44 !! * Substitutions 38 45 # include "domzgr_substitute.h90" 46 # include "vectopt_loop_substitute.h90" 39 47 !!---------------------------------------------------------------------- 40 48 !! NEMO/OPA 4.0 , NEMO Consortium (2011) … … 51 59 !! 52 60 !! ** Method : - Read namelist namsbc_wave 53 !! - Read Cd_n10 fields in netcdf files 54 !! - Read stokes drift 2d in netcdf files 55 !! - Read wave number in netcdf files 56 !! - Compute 3d stokes drift using monochromatic 57 !! ** action : 58 !! 59 !!--------------------------------------------------------------------- 60 USE oce, ONLY : un,vn,hdivn,rotn 61 USE divcur 62 USE wrk_nemo 63 #if defined key_bdy 64 USE bdy_oce, ONLY : bdytmask 65 #endif 66 INTEGER, INTENT( in ) :: kt ! ocean time step 61 !! - Read Cd_n10 fields in netcdf files 62 !! ** action : 63 !! 64 !!--------------------------------------------------------------------- 65 INTEGER, INTENT( in ) :: kt ! ocean time step 67 66 INTEGER :: ierror ! return error code 68 INTEGER :: ifpr, jj,ji,jk 69 INTEGER :: ios ! Local integer output status for namelist read 70 REAL(wp),DIMENSION(:,:,:),POINTER :: udummy,vdummy,hdivdummy,rotdummy 71 REAL :: z2dt,z1_2dt 72 TYPE(FLD_N), DIMENSION(jpfld) :: slf_i ! array of namelist informations on the fields to read 73 CHARACTER(len=100) :: cn_dir ! Root directory for location of drag coefficient files 74 TYPE(FLD_N) :: sn_cdg, sn_usd, sn_vsd, sn_wn ! informations about the fields to be read 75 !!--------------------------------------------------------------------- 76 NAMELIST/namsbc_wave/ sn_cdg, cn_dir, sn_usd, sn_vsd, sn_wn 67 CHARACTER(len=100) :: cn_dir_cdg ! Rootdirectory for location of drag coefficient files 68 TYPE(FLD_N) :: sn_cdg ! informationsabout the fields to be read 69 !!--------------------------------------------------------------------- 70 NAMELIST/namsbc_wave/ sn_cdg, cn_dir_cdg 77 71 !!--------------------------------------------------------------------- 78 72 … … 83 77 IF( kt == nit000 ) THEN ! First call kt=nit000 ! 84 78 ! ! -------------------- ! 85 REWIND( numnam_ref ) ! Namelist namsbc_wave in reference namelist : File for drag coeff. from wave model 86 READ ( numnam_ref, namsbc_wave, IOSTAT = ios, ERR = 901) 87 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in reference namelist', lwp ) 88 89 REWIND( numnam_cfg ) ! Namelist namsbc_wave in configuration namelist : File for drag coeff. from wave model 90 READ ( numnam_cfg, namsbc_wave, IOSTAT = ios, ERR = 902 ) 91 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in configuration namelist', lwp ) 92 IF(lwm) WRITE ( numond, namsbc_wave ) 93 ! 94 95 IF ( ln_cdgw ) THEN 96 ALLOCATE( sf_cd(1), STAT=ierror ) !* allocate and fill sf_wave with sn_cdg 97 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 98 ! 99 ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1) ) 100 IF( sn_cdg%ln_tint ) ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) ) 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.0 104 ENDIF 105 IF ( ln_sdw ) THEN 106 slf_i(jp_usd) = sn_usd ; slf_i(jp_vsd) = sn_vsd; slf_i(jp_wn) = sn_wn 107 ALLOCATE( sf_sd(3), STAT=ierror ) !* allocate and fill sf_wave with sn_cdg 108 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 109 ! 110 DO ifpr= 1, jpfld 111 ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) ) 112 IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) ) 113 END DO 114 CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 115 ALLOCATE( usd2d(jpi,jpj),vsd2d(jpi,jpj),uwavenum(jpi,jpj),vwavenum(jpi,jpj) ) 116 ALLOCATE( usd3d(jpi,jpj,jpk),vsd3d(jpi,jpj,jpk),wsd3d(jpi,jpj,jpk) ) 117 usd2d(:,:) = 0.0 ; vsd2d(:,:) = 0.0 ; uwavenum(:,:) = 0.0 ; vwavenum(:,:) = 0.0 118 usd3d(:,:,:) = 0.0 ;vsd3d(:,:,:) = 0.0 ; wsd3d(:,:,:) = 0.0 79 ! !* set file information 80 ! (default values) 81 ! ... default values (NB: frequency positive => hours, negative => 82 ! months) 83 ! ! file ! frequency ! variable ! time intep ! clim 84 ! ! 'yearly' or ! weights ! rotation ! 85 ! ! name ! (hours) ! name ! (T/F) ! 86 ! (T/F) ! 'monthly' ! filename ! pairs ! 87 sn_cdg = FLD_N('cdg_wave' , 1 ,'drag_coeff', .true. ,.false. , 'daily' , '' , '','' ) 88 cn_dir_cdg = './' ! directory in which the Patm data are 89 90 91 REWIND( numnam_ref ) !* read in namlist namsbc_wave 92 READ ( numnam_ref, namsbc_wave ) 93 REWIND( numnam_ref ) 94 REWIND( numnam_cfg ) !* read in namlist namsbc_wave 95 READ ( numnam_cfg, namsbc_wave ) 96 REWIND( numnam_cfg ) 97 ! 98 99 ALLOCATE( sf_wave(1), STAT=ierror ) !* allocate and fillsf_wave with sn_cdg 100 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocatesf_wave structure' ) 101 ! 102 CALL fld_fill( sf_wave, (/ sn_cdg /), cn_dir_cdg, 'sbc_wave', 'Wavemodule ', 'namsbc_wave' ) 103 ALLOCATE( sf_wave(1)%fnow(jpi,jpj,1) ) 104 IF( sn_cdg%ln_tint ) ALLOCATE( sf_wave(1)%fdta(jpi,jpj,1,2) ) 105 !ALLOCATE( cdn_wave(jpi,jpj) ) 106 ! Allocation done by sbc_oce 107 cdn_wave(:,:) = 0.0 108 ENDIF 109 ! 110 ! 111 CALL fld_read( kt, nn_fsbc, sf_wave ) !* read drag coefficient fromexternal forcing 112 cdn_wave(:,:) = sf_wave(1)%fnow(:,:,1) 113 END SUBROUTINE sbc_wave 114 115 116 SUBROUTINE sbc_wavepar( kt ) 117 !!--------------------------------------------------------------------- 118 !! *** ROUTINE sbc_wavepar *** 119 !! 120 !! ** Purpose : Provide at each time step wave model parameters, including the 121 !! Stokes drift (east and north components), significant wave height and the 122 !! mean wave period as well as the normalized stress and energy flux into the 123 !! ocean for TKE. 124 !! 125 !! 126 !! ** Method : - Read namelist namsbc_wavepar 127 !! - Read fields in NetCDF files 128 !! ** action : 129 !! 130 !!--------------------------------------------------------------------- 131 INTEGER, INTENT( in ) :: kt ! ocean time step 132 !! 133 INTEGER :: ierror ! return error code 134 INTEGER :: ji,iii,jjj ! dummy loop index, Victor added iii, jjj 135 INTEGER :: jfld ! dummy loop arguments 136 !! 137 CHARACTER(len=100) :: cn_dir_wavepar ! Root directory for location of ECWAM wave parameter fields 138 TYPE(FLD_N), DIMENSION(jpfld_wavepar) :: slf_i ! array of namelist informations on the fields to read 139 TYPE(FLD_N) :: sn_ust, sn_vst, sn_swh, sn_mwp ! information about the fields to be read 140 TYPE(FLD_N) :: sn_phioc, sn_tauoc 141 TYPE(FLD_N) :: sn_wspd 142 !TYPE(FLD_N) :: sn_cdww 143 !!--------------------------------------------------------------------- 144 !NAMELIST/namsbc_wavepar/ sn_ust, sn_vst, sn_swh, sn_mwp, sn_phioc, sn_tauoc, cn_dir_wavepar 145 NAMELIST/namsbc_wavepar/ sn_ust, sn_vst, sn_swh, sn_mwp, sn_wspd, sn_phioc, sn_tauoc, cn_dir_wavepar 146 !!--------------------------------------------------------------------- 147 148 !!---------------------------------------------------------------------- 149 ! 150 ! 151 ! ! -------------------- ! 152 IF( kt == nit000 ) THEN ! First call kt=nit000 ! 153 ! ! -------------------- ! 154 ! set file information (default values) 155 ! ... default values (NB: frequency positive => hours, negative => months) 156 ! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! weights ! rotation ! 157 ! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! filename ! pairs ! 158 sn_ust = FLD_N( 'ust' , 6 , 'ust' , .true. , .false. , 'monthly' , '' , '' ,'' ) 159 sn_vst = FLD_N( 'vst' , 6 , 'vst' , .true. , .false. , 'monthly' , '' , '' ,'' ) 160 sn_swh = FLD_N( 'swh' , 6 , 'swh' , .true. , .false. , 'monthly' , '' , '' ,'' ) 161 sn_mwp = FLD_N( 'mwp' , 6 , 'mwp' , .true. , .false. , 'monthly' , '' , '' ,'' ) 162 !sn_cdww = FLD_N( 'cdww' , 6 , 'cdww' , .true. , .false. , 'monthly' , '' , '' ) 163 sn_wspd = FLD_N( 'wspd' , 6 , 'wspd' , .true. , .false. , 'monthly' , '' , '','' ) 164 sn_phioc = FLD_N( 'phioc' , 6 , 'phioc' , .true. , .false. , 'monthly' , '' , '','' ) 165 sn_tauoc = FLD_N( 'tauoc' , 6 , 'tauoc' , .true. , .false. , 'monthly' , '' , '','' ) 166 cn_dir_wavepar = './' ! directory in which the wave data are found 167 168 REWIND( numnam_ref ) !* read in namlist namsbc_wave 169 READ ( numnam_ref, namsbc_wavepar ) 170 REWIND ( numnam_ref ) 171 172 REWIND( numnam_cfg ) !* read in namlist namsbc_wave 173 READ ( numnam_cfg, namsbc_wavepar ) 174 REWIND ( numnam_cfg ) 175 ! 176 177 178 slf_i(jp_ust) = sn_ust 179 slf_i(jp_vst) = sn_vst 180 slf_i(jp_swh) = sn_swh 181 slf_i(jp_mwp) = sn_mwp 182 !slf_i(jp_cdww) = sn_cdww 183 slf_i(jp_wspd) = sn_wspd 184 slf_i(jp_phioc) = sn_phioc 185 slf_i(jp_tauoc) = sn_tauoc 186 187 ALLOCATE( sf_wavepar(jpfld_wavepar), STAT=ierror ) !* set sf_wavepar structure 188 IF ( ierror > 0 ) THEN 189 CALL ctl_stop( 'STOP', 'sbc_wavepar: unable to allocate sf_wavepar structure' ) ; RETURN 119 190 ENDIF 191 ! 192 jfld = jpfld_wavepar 193 DO ji = 1, jpfld_wavepar 194 ALLOCATE( sf_wavepar(ji)%fnow(jpi,jpj,1) ) 195 IF ( slf_i(ji)%ln_tint ) ALLOCATE( sf_wavepar(ji)%fdta(jpi,jpj,1,2) ) 196 ENDDO 197 ALLOCATE( ust_wavepar(jpi,jpj) ) 198 ALLOCATE( vst_wavepar(jpi,jpj) ) 199 ALLOCATE( swh_wavepar(jpi,jpj) ) 200 ALLOCATE( mwp_wavepar(jpi,jpj) ) 201 !ALLOCATE( cdww_wavepar(jpi,jpj) ) 202 !ALLOCATE( wspd_wavepar(jpi,jpj) ) 203 ALLOCATE( phioc_wavepar(jpi,jpj) ) 204 ALLOCATE( tauoc_wavepar(jpi,jpj) ) 205 ! 206 CALL fld_fill( sf_wavepar, slf_i, cn_dir_wavepar, 'sbc_wavepar', 'Wave module ', 'namsbc_wavepar' ) 207 ! 120 208 ENDIF 121 209 ! 122 ! 123 IF ( ln_cdgw ) THEN 124 CALL fld_read( kt, nn_fsbc, sf_cd ) !* read drag coefficient from external forcing 125 cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1) 126 ENDIF 127 IF ( ln_sdw ) THEN 128 CALL fld_read( kt, nn_fsbc, sf_sd ) !* read drag coefficient from external forcing 129 130 ! Interpolate wavenumber, stokes drift into the grid_V and grid_V 131 !------------------------------------------------- 132 133 DO jj = 1, jpjm1 134 DO ji = 1, jpim1 135 uwavenum(ji,jj)=0.5 * ( 2. - umask(ji,jj,1) ) * ( sf_sd(3)%fnow(ji,jj,1) * tmask(ji,jj,1) & 136 & + sf_sd(3)%fnow(ji+1,jj,1) * tmask(ji+1,jj,1) ) 137 138 vwavenum(ji,jj)=0.5 * ( 2. - vmask(ji,jj,1) ) * ( sf_sd(3)%fnow(ji,jj,1) * tmask(ji,jj,1) & 139 & + sf_sd(3)%fnow(ji,jj+1,1) * tmask(ji,jj+1,1) ) 140 141 usd2d(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( sf_sd(1)%fnow(ji,jj,1) * tmask(ji,jj,1) & 142 & + sf_sd(1)%fnow(ji+1,jj,1) * tmask(ji+1,jj,1) ) 143 144 vsd2d(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( sf_sd(2)%fnow(ji,jj,1) * tmask(ji,jj,1) & 145 & + sf_sd(2)%fnow(ji,jj+1,1) * tmask(ji,jj+1,1) ) 146 END DO 147 END DO 148 149 !Computation of the 3d Stokes Drift 150 DO jk = 1, jpk 151 DO jj = 1, jpj-1 152 DO ji = 1, jpi-1 153 usd3d(ji,jj,jk) = usd2d(ji,jj)*exp(2.0*uwavenum(ji,jj)*(-MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj ,jk)))) 154 vsd3d(ji,jj,jk) = vsd2d(ji,jj)*exp(2.0*vwavenum(ji,jj)*(-MIN( gdept_0(ji,jj,jk) , gdept_0(ji ,jj+1,jk)))) 155 END DO 156 END DO 157 usd3d(jpi,:,jk) = usd2d(jpi,:)*exp( 2.0*uwavenum(jpi,:)*(-gdept_0(jpi,:,jk)) ) 158 vsd3d(:,jpj,jk) = vsd2d(:,jpj)*exp( 2.0*vwavenum(:,jpj)*(-gdept_0(:,jpj,jk)) ) 159 END DO 160 161 CALL wrk_alloc( jpi,jpj,jpk,udummy,vdummy,hdivdummy,rotdummy) 162 163 udummy(:,:,:)=un(:,:,:) 164 vdummy(:,:,:)=vn(:,:,:) 165 hdivdummy(:,:,:)=hdivn(:,:,:) 166 rotdummy(:,:,:)=rotn(:,:,:) 167 un(:,:,:)=usd3d(:,:,:) 168 vn(:,:,:)=vsd3d(:,:,:) 169 CALL div_cur(kt) 170 ! !------------------------------! 171 ! ! Now Vertical Velocity ! 172 ! !------------------------------! 173 z2dt = 2._wp * rdt ! set time step size (Euler/Leapfrog) 174 175 z1_2dt = 1.e0 / z2dt 176 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 177 ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise 178 wsd3d(:,:,jk) = wsd3d(:,:,jk+1) - fse3t_n(:,:,jk) * hdivn(:,:,jk) & 179 & - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) & 180 & * tmask(:,:,jk) * z1_2dt 181 #if defined key_bdy 182 wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:) 183 #endif 184 END DO 185 hdivn(:,:,:)=hdivdummy(:,:,:) 186 rotn(:,:,:)=rotdummy(:,:,:) 187 vn(:,:,:)=vdummy(:,:,:) 188 un(:,:,:)=udummy(:,:,:) 189 CALL wrk_dealloc( jpi,jpj,jpk,udummy,vdummy,hdivdummy,rotdummy) 190 ENDIF 191 END SUBROUTINE sbc_wave 210 CALL fld_read( kt, nn_fsbc, sf_wavepar ) !* read wave parameters from ECWAM NetCDF file 211 ust_wavepar(:,:) = sf_wavepar(jp_ust)%fnow(:,:,1) 212 vst_wavepar(:,:) = sf_wavepar(jp_vst)%fnow(:,:,1) 213 swh_wavepar(:,:) = sf_wavepar(jp_swh)%fnow(:,:,1) 214 mwp_wavepar(:,:) = sf_wavepar(jp_mwp)%fnow(:,:,1) 215 !cdww_wavepar(:,:) = sf_wavepar(jp_cdww)%fnow(:,:,1) 216 wspd_wavepar(:,:) = sf_wavepar(jp_wspd)%fnow(:,:,1) 217 phioc_wavepar(:,:) = sf_wavepar(jp_phioc)%fnow(:,:,1) 218 tauoc_wavepar(:,:) = sf_wavepar(jp_tauoc)%fnow(:,:,1) 219 220 221 rn_crban(:,:)=29*phioc_wavepar(:,:) ! Alfa is phioc*sqrt(rw/ra)sbc_wa 222 ! Limit cr_ban between 10 and 300 223 DO iii=1,jpi 224 DO jjj=1,jpj 225 226 IF (rn_crban(iii,jjj) < 10 ) THEN 227 rn_crban(iii,jjj)=10 228 ELSEIF (rn_crban(iii,jjj) > 300) THEN 229 rn_crban(iii,jjj)=300 230 ENDIF 231 232 ENDDO 233 ENDDO 234 235 END SUBROUTINE sbc_wavepar 192 236 193 !!====================================================================== 237 194 238 END MODULE sbcwave -
branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
r6204 r7807 24 24 USE phycst ! physical constants 25 25 USE zdfmxl ! mixed layer 26 USE sbcwave 27 ! 26 28 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 27 29 USE lib_mpp ! MPP manager … … 47 49 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustarb2 !: Squared bottom velocity scale at T-points 48 50 51 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rsbc_tke1 52 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rsbc_tke3 53 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rsbc_psi1 54 49 55 ! !! ** Namelist namzdf_gls ** 56 LOGICAL :: ln_crban = .FALSE. !! WAVE2NEMO added ln_crban (=.TRUE. use Craig and Banner scheme) 50 57 LOGICAL :: ln_length_lim ! use limit on the dissipation rate under stable stratification (Galperin et al. 1988) 51 58 LOGICAL :: ln_sigpsi ! Activate Burchard (2003) modification for k-eps closure & wave breaking mixing … … 59 66 REAL(wp) :: rn_emin ! minimum value of TKE (m2/s2) 60 67 REAL(wp) :: rn_charn ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value) 61 REAL(wp) :: rn_crban 68 REAL(wp) :: rn_crban_default ! Craig and Banner constant for surface breaking waves mixing 62 69 REAL(wp) :: rn_hsro ! Minimum surface roughness 63 REAL(wp) :: rn_frac_hs ! Fraction of wave height as surface roughness (if nn_z0_met >1)70 REAL(wp) :: rn_frac_hs ! Fraction of wave height as surface roughness (if nn_z0_met = 1) 64 71 65 72 REAL(wp) :: rcm_sf = 0.73_wp ! Shear free turbulence parameters … … 92 99 REAL(wp) :: rtrans = 0.1_wp 93 100 REAL(wp) :: rc02, rc02r, rc03, rc04 ! coefficients deduced from above parameters 94 REAL(wp) :: rsbc_tke1 , rsbc_tke2, rfact_tke! - - - -95 REAL(wp) :: rsbc_psi 1, rsbc_psi2, rfact_psi ! - - - -101 REAL(wp) :: rsbc_tke1_default, rsbc_tke2, rfact_tke ! - - - - 102 REAL(wp) :: rsbc_psi2, rsbc_psi3, rfact_psi ! - - - - 96 103 REAL(wp) :: rsbc_zs1, rsbc_zs2 ! - - - - 97 104 REAL(wp) :: rc0, rc2, rc3, rf6, rcff, rc_diff ! - - - - … … 115 122 !! *** FUNCTION zdf_gls_alloc *** 116 123 !!---------------------------------------------------------------------- 117 ALLOCATE( mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) , & 118 & ustars2(jpi,jpj) , ustarb2(jpi,jpj) , STAT= zdf_gls_alloc ) 124 ALLOCATE( mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) , & 125 & ustars2(jpi,jpj), ustarb2(jpi,jpj),rsbc_tke1(jpi,jpj), & 126 & rsbc_tke3(jpi,jpj), rsbc_psi1(jpi,jpj) , STAT=zdf_gls_alloc ) 119 127 ! 120 128 IF( lk_mpp ) CALL mpp_sum ( zdf_gls_alloc ) … … 161 169 ustars2 = 0._wp ; ustarb2 = 0._wp ; psi = 0._wp ; zwall_psi = 0._wp 162 170 171 172 ! variable initialization 173 IF( ln_crban ) THEN 174 rsbc_tke1(:,:) = (-rsc_tke*rn_crban(:,:)/(rcm_sf*ra_sf*rl_sf))**(2._wp/3._wp) ! k_eps = 53.Dirichlet + Wave breaking 175 rsbc_tke3(:,:) = rdt * rn_crban(:,:) ! Neumann + Wave breaking 176 rsbc_psi1(:,:) = rc0**rpp * rsbc_tke1(:,:)**rmm * rl_sf**rnn ! Dirichlet + Wave breaking 177 ENDIF 178 163 179 IF( kt /= nit000 ) THEN ! restore before value to compute tke 164 180 avt (:,:,:) = avt_k (:,:,:) … … 197 213 zdep(:,:) = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(:,:),rsmall)))) ! Wave age (eq. 10) 198 214 zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 199 ! 215 CASE ( 3 ) ! Roughness given by the wave model (coupled or read in file) 216 WHERE( swh_wavepar == 0._wp ) ! surface roughness length according to Charnock formula when sign. wave height 0 217 zhsro = MAX(rn_charn / grav * ustars2, rn_hsro) 218 ELSEWHERE 219 zhsro = MAX(swh_wavepar, rn_hsro) 220 END WHERE 200 221 END SELECT 201 222 … … 314 335 ! 315 336 CASE ( 0 ) ! Dirichlet case 316 ! First level 317 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 318 en(:,:,1) = MAX(en(:,:,1), rn_emin) 319 z_elem_a(:,:,1) = en(:,:,1) 320 z_elem_c(:,:,1) = 0._wp 321 z_elem_b(:,:,1) = 1._wp 322 ! 323 ! One level below 324 en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2)) & 325 & / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 326 en(:,:,2) = MAX(en(:,:,2), rn_emin ) 327 z_elem_a(:,:,2) = 0._wp 328 z_elem_c(:,:,2) = 0._wp 329 z_elem_b(:,:,2) = 1._wp 330 ! 331 ! 337 IF( ln_crban ) THEN ! wave induced mixing case with forced/coupled fields 338 ! First level 339 en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin ) 340 z_elem_a(:,:,1) = en(:,:,1) 341 z_elem_c(:,:,1) = 0._wp 342 z_elem_b(:,:,1) = 1._wp 343 344 ! One level below 345 en(:,:,2) = MAX( rsbc_tke1(:,:) * ustars2(:,:) * ((zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin ) 346 z_elem_a(:,:,2) = 0._wp 347 z_elem_c(:,:,2) = 0._wp 348 z_elem_b(:,:,2) = 1._wp 349 ELSE ! wave induced mixing case with default values 350 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default)**(2._wp/3._wp) 351 en(:,:,1) = MAX(en(:,:,1), rn_emin) 352 z_elem_a(:,:,1) = en(:,:,1) 353 z_elem_c(:,:,1) = 0._wp 354 z_elem_b(:,:,1) = 1._wp 355 ! 356 ! One level below 357 en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default * ((zhsro(:,:)+fsdepw(:,:,2)) & 358 & / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 359 en(:,:,2) = MAX(en(:,:,2), rn_emin ) 360 z_elem_a(:,:,2) = 0._wp 361 z_elem_c(:,:,2) = 0._wp 362 z_elem_b(:,:,2) = 1._wp 363 ! 364 ! 365 ENDIF 332 366 CASE ( 1 ) ! Neumann boundary condition on d(e)/dz 333 ! 334 ! Dirichlet conditions at k=1 335 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 336 en(:,:,1) = MAX(en(:,:,1), rn_emin) 337 z_elem_a(:,:,1) = en(:,:,1) 338 z_elem_c(:,:,1) = 0._wp 339 z_elem_b(:,:,1) = 1._wp 340 ! 341 ! at k=2, set de/dz=Fw 342 !cbr 343 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 344 z_elem_a(:,:,2) = 0._wp 345 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:)) )) 346 zflxs(:,:) = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) & 347 & * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf) 348 349 en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2) 350 ! 351 ! 367 IF( ln_crban ) THEN ! Shear free case: d(e)/dz=Fw with forced/coupled fields 368 ! Dirichlet conditions at k=1 369 en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin ) 370 z_elem_a(:,:,1) = en(:,:,1) 371 z_elem_c(:,:,1) = 0._wp 372 z_elem_b(:,:,1) = 1._wp 373 ! at k=2, set de/dz=Fw 374 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 375 z_elem_a(:,:,2) = 0._wp 376 zflxs(:,:) = rsbc_tke3(:,:) * ustars2(:,:)**1.5_wp * ((zhsro(:,:)+fsdept(:,:,1) ) / zhsro(:,:) )**(1.5*ra_sf) 377 en(:,:,2) = en(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 378 ELSE ! Shear free case: d(e)/dz=Fw with default values 379 ! Dirichlet conditions at k=1 380 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default)**(2._wp/3._wp) 381 en(:,:,1) = MAX(en(:,:,1), rn_emin) 382 z_elem_a(:,:,1) = en(:,:,1) 383 z_elem_c(:,:,1) = 0._wp 384 z_elem_b(:,:,1) = 1._wp 385 ! 386 ! at k=2, set de/dz=Fw 387 !cbr 388 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 389 z_elem_a(:,:,2) = 0._wp 390 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:)) )) 391 zflxs(:,:) = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) & 392 & * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf) 393 394 en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2) 395 ! 396 ! 397 ENDIF 352 398 END SELECT 353 399 … … 536 582 ! 537 583 CASE ( 0 ) ! Dirichlet boundary conditions 538 ! 539 ! Surface value 540 zdep(:,:) = zhsro(:,:) * rl_sf ! Cosmetic 541 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 542 z_elem_a(:,:,1) = psi(:,:,1) 543 z_elem_c(:,:,1) = 0._wp 544 z_elem_b(:,:,1) = 1._wp 545 ! 546 ! One level below 547 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdepw(:,:,2)/zhsro(:,:) ))) 548 zdep(:,:) = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:) 549 psi (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 550 z_elem_a(:,:,2) = 0._wp 551 z_elem_c(:,:,2) = 0._wp 552 z_elem_b(:,:,2) = 1._wp 553 ! 554 ! 584 IF( ln_crban ) THEN ! Wave induced mixing case 585 ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2 586 ! balance between the production and the 587 ! dissipation terms including the wave effect 588 ! Surface value 589 zdep(:,:) = zhsro(:,:) * rl_sf ! Cosmetic 590 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 591 z_elem_a(:,:,1) = psi(:,:,1) 592 z_elem_c(:,:,1) = 0._wp 593 z_elem_b(:,:,1) = 1._wp 594 ! 595 ! One level below 596 zex1 = (rmm*ra_sf+rnn) 597 zex2 = (rmm*ra_sf) 598 zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**zex1 ) / zhsro(:,:)**zex2 599 psi (:,:,2) = rsbc_psi1(:,:) * ustars2(:,:)**rmm * zdep(:,:) * tmask(:,:,1) 600 z_elem_a(:,:,2) = 0._wp 601 z_elem_c(:,:,2) = 0._wp 602 z_elem_b(:,:,2) = 1._wp 603 ! 604 ! 605 ELSE ! Wave induced mixing case with default values 606 ! Surface value 607 zdep(:,:) = zhsro(:,:) * rl_sf ! Cosmetic 608 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 609 z_elem_a(:,:,1) = psi(:,:,1) 610 z_elem_c(:,:,1) = 0._wp 611 z_elem_b(:,:,1) = 1._wp 612 ! 613 ! One level below 614 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdepw(:,:,2)/zhsro(:,:) ))) 615 zdep(:,:) = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:) 616 psi (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 617 z_elem_a(:,:,2) = 0._wp 618 z_elem_c(:,:,2) = 0._wp 619 z_elem_b(:,:,2) = 1._wp 620 ! 621 ! 622 ENDIF 555 623 CASE ( 1 ) ! Neumann boundary condition on d(psi)/dz 556 ! 557 ! Surface value: Dirichlet 558 zdep(:,:) = zhsro(:,:) * rl_sf 559 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 560 z_elem_a(:,:,1) = psi(:,:,1) 561 z_elem_c(:,:,1) = 0._wp 562 z_elem_b(:,:,1) = 1._wp 563 ! 564 ! Neumann condition at k=2 565 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 566 z_elem_a(:,:,2) = 0._wp 567 ! 568 ! Set psi vertical flux at the surface: 569 zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 570 zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 571 zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 572 zdep(:,:) = rsbc_psi1 * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 573 & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.) 574 zflxs(:,:) = zdep(:,:) * zflxs(:,:) 575 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 576 577 ! 578 ! 624 IF( ln_crban ) THEN ! Wave induced mixing case with forced/coupled fields 625 ! 626 zdep(:,:) = rl_sf * zhsro(:,:) 627 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 628 z_elem_a(:,:,1) = psi(:,:,1) 629 z_elem_c(:,:,1) = 0._wp 630 z_elem_b(:,:,1) = 1._wp 631 ! 632 ! Neumann condition at k=2 633 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 634 z_elem_a(:,:,2) = 0._wp 635 ! 636 ! Set psi vertical flux at the surface: 637 zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(rmm*ra_sf+rnn-1._wp) / zhsro(:,:)**(rmm*ra_sf) 638 zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) & 639 & * en(:,:,1)**rmm * zdep 640 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 641 ! 642 ELSE ! Wave induced mixing case with default values 643 ! Surface value: Dirichlet 644 zdep(:,:) = zhsro(:,:) * rl_sf 645 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 646 z_elem_a(:,:,1) = psi(:,:,1) 647 z_elem_c(:,:,1) = 0._wp 648 z_elem_b(:,:,1) = 1._wp 649 ! 650 ! Neumann condition at k=2 651 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 652 z_elem_a(:,:,2) = 0._wp 653 ! 654 ! Set psi vertical flux at the surface: 655 zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 656 zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 657 zflxs(:,:) = (rnn + rsbc_tke1_default * (rnn + rmm*ra_sf) * zdep(:,:)) * & 658 (1._wp + rsbc_tke1_default * zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 659 zdep(:,:) = rsbc_psi1(:,:) * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 660 & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.) 661 zflxs(:,:) = zdep(:,:) * zflxs(:,:) 662 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 663 ! 664 ! 665 ENDIF 579 666 END SELECT 580 667 … … 866 953 !! 867 954 NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, & 868 & rn_clim_galp, ln_ sigpsi, rn_hsro, &869 & rn_crban , rn_charn, rn_frac_hs,&955 & rn_clim_galp, ln_crban, ln_sigpsi, rn_hsro, & 956 & rn_crban_default, rn_charn, rn_frac_hs,& 870 957 & nn_bc_surf, nn_bc_bot, nn_z0_met, & 871 958 & nn_stab_func, nn_clos … … 895 982 WRITE(numout,*) ' TKE Bottom boundary condition nn_bc_bot = ', nn_bc_bot 896 983 WRITE(numout,*) ' Modify psi Schmidt number (wb case) ln_sigpsi = ', ln_sigpsi 897 WRITE(numout,*) ' Craig and Banner coefficient rn_crban = ', rn_crban 984 WRITE(numout,*) ' Craig and Banner coefficient (default) rn_crban = ', rn_crban_default 985 WRITE(numout,*) ' Craig and Banner scheme ln_crban = ', ln_crban 898 986 WRITE(numout,*) ' Charnock coefficient rn_charn = ', rn_charn 899 987 WRITE(numout,*) ' Surface roughness formula nn_z0_met = ', nn_z0_met … … 909 997 910 998 ! !* Check of some namelist values 911 IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' ) 912 IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' ) 913 IF( nn_z0_met < 0 .OR. nn_z0_met > 2 ) CALL ctl_stop( 'bad flag: nn_z0_met is 0, 1 or 2' ) 914 IF( nn_stab_func < 0 .OR. nn_stab_func > 3 ) CALL ctl_stop( 'bad flag: nn_stab_func is 0, 1, 2 and 3' ) 915 IF( nn_clos < 0 .OR. nn_clos > 3 ) CALL ctl_stop( 'bad flag: nn_clos is 0, 1, 2 or 3' ) 999 IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' ) 1000 IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' ) 1001 IF( nn_z0_met < 0 .OR. nn_z0_met > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_z0_met is 0, 1, 2 or 3' ) 1002 ! IF( nn_z0_met == 3 .AND. .NOT.ln_sdw ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_sdw=T' ) 1003 IF( nn_stab_func < 0 .OR. nn_stab_func > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_stab_func is 0, 1, 2 and 3' ) 1004 IF( nn_clos < 0 .OR. nn_clos > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' ) 916 1005 917 1006 SELECT CASE ( nn_clos ) !* set the parameters for the chosen closure … … 1085 1174 & - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) ) 1086 1175 1087 IF ( rn_crban==0._wp ) THEN1176 IF( .NOT. ln_crban .AND. rn_crban_default==0._wp ) THEN 1088 1177 rl_sf = vkarmn 1089 1178 ELSE … … 1124 1213 rc03 = rc02 * rc0 1125 1214 rc04 = rc03 * rc0 1126 rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf ! Dirichlet + Wave breaking1127 rsbc_tke2 = rdt * rn_crban / rl_sf ! Neumann + Wave breaking1128 zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp )1129 rtrans = 0.2_wp / zcr ! Ad. inverse transition length between log and wave layer1215 rsbc_tke1_default = -3._wp/2._wp*rn_crban_default*ra_sf*rl_sf 1216 rsbc_tke2 = rdt * rn_crban_default / rl_sf 1217 zcr = MAX( rsmall, rsbc_tke1_default**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 1218 rtrans = 0.2_wp / zcr 1130 1219 rsbc_zs1 = rn_charn/grav ! Charnock formula for surface roughness 1131 1220 rsbc_zs2 = rn_frac_hs / 0.85_wp / grav * 665._wp ! Rascle formula for surface roughness -
branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/step.F90
r6204 r7807 323 323 CALL dyn_vor( kstp ) ! vorticity term including Coriolis 324 324 CALL dyn_ldf( kstp ) ! lateral mixing 325 IF ( ln_stcor ) CALL dyn_stcor ( kstp ) ! WAVE2NEMO add Stokes-Coriolis forcing 325 326 IF( ln_neptsimp ) CALL dyn_nept_cor( kstp ) ! add Neptune velocities (simplified) 326 327 #if defined key_agrif -
branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
r6204 r7807 110 110 USE timing ! Timing 111 111 112 USE dynstcor ! WAVE2NEMO Stokes-Coriolis forcing 112 113 #if defined key_agrif 113 114 USE agrif_opa_sponge ! Momemtum and tracers sponges
Note: See TracChangeset
for help on using the changeset viewer.
