- Timestamp:
- 2017-12-13T15:58:53+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
r7864 r9019 19 19 USE oce ! ocean variables 20 20 USE sbc_oce ! Surface boundary condition: ocean fields 21 USE zdf_oce, ONLY : ln_zdf qiao21 USE zdf_oce, ONLY : ln_zdfswm 22 22 USE bdy_oce ! open boundary condition variables 23 23 USE domvvl ! domain: variable volume layers … … 87 87 INTEGER :: jj, ji, jk ! dummy loop argument 88 88 INTEGER :: ik ! local integer 89 REAL(wp) :: ztransp, zfac, ztemp, zsp0 90 REAL(wp) :: zdep_u, zdep_v, zkh_u, zkh_v, zda_u, zda_v 89 REAL(wp) :: ztransp, zfac, zsp0 90 REAL(wp) :: zdepth, zsqrt_depth, zexp_depth, z_two_thirds, zsqrtpi !sqrt of pi 91 REAL(wp) :: zbot_u, zbot_v, zkb_u, zkb_v, zke3_u, zke3_v, zda_u, zda_v 92 REAL(wp) :: zstokes_psi_u_bot, zstokes_psi_v_bot 91 93 REAL(wp), DIMENSION(:,:) , POINTER :: zk_t, zk_u, zk_v, zu0_sd, zv0_sd ! 2D workspace 94 REAL(wp), DIMENSION(:,:) , POINTER :: zstokes_psi_u_top, zstokes_psi_v_top ! 2D workspace 92 95 REAL(wp), DIMENSION(:,:,:), POINTER :: ze3divh ! 3D workspace 93 96 !!--------------------------------------------------------------------- … … 95 98 CALL wrk_alloc( jpi,jpj,jpk, ze3divh ) 96 99 CALL wrk_alloc( jpi,jpj, zk_t, zk_u, zk_v, zu0_sd, zv0_sd ) 97 ! 98 ! 100 CALL wrk_alloc( jpi,jpj, zstokes_psi_u_top, zstokes_psi_v_top) 101 ! 102 ! 103 zsqrtpi = SQRT(rpi) 104 z_two_thirds = 2.0_wp / 3.0_wp 99 105 zfac = 2.0_wp * rpi / 16.0_wp 100 106 DO jj = 1, jpj ! exp. wave number at t-point (Eq. (19) in Breivick et al. (2014) ) … … 106 112 tsd2d(ji,jj) = zsp0 107 113 ! Wavenumber scale 108 zk_t(ji,jj) = ABS( zsp0 ) / MAX( ABS( 5.97_wp*ztransp ) , 0.0000001_wp ) 109 END DO 110 END DO 114 zk_t(ji,jj) = (1.0_wp-2.0_wp/3.0_wp)*zsp0/MAX(2.0_wp*ztransp,0.0000001_wp) 115 END DO 116 END DO 117 ! 111 118 DO jj = 1, jpjm1 ! exp. wave number & Stokes drift velocity at u- & v-points 112 119 DO ji = 1, jpim1 … … 120 127 ! 121 128 ! !== horizontal Stokes Drift 3D velocity ==! 129 DO jj = 1, jpjm1 ! exp. wave number & Stokes drift velocity at u- & v-points 130 DO ji = 1, jpim1 131 zstokes_psi_u_top(ji,jj) = 0._wp 132 zstokes_psi_v_top(ji,jj) = 0._wp 133 END DO 134 END DO 135 122 136 DO jk = 1, jpkm1 123 137 DO jj = 2, jpjm1 124 138 DO ji = 2, jpim1 125 zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) ) 126 zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) ) 127 ! 128 zkh_u = zk_u(ji,jj) * zdep_u ! k * depth 129 zkh_v = zk_v(ji,jj) * zdep_v 130 ! ! Depth attenuation 131 zda_u = EXP( -2.0_wp*zkh_u ) / ( 1.0_wp + 8.0_wp*zkh_u ) 132 zda_v = EXP( -2.0_wp*zkh_v ) / ( 1.0_wp + 8.0_wp*zkh_v ) 139 zbot_u = 0.5_wp * ( gdepw_n(ji,jj,jk+1) + gdepw_n(ji+1,jj,jk+1) ) 140 zbot_v = 0.5_wp * ( gdepw_n(ji,jj,jk+1) + gdepw_n(ji,jj+1,jk+1) ) 141 zkb_u = 2.0_wp * zk_u(ji,jj) * zbot_u ! 2k * bottom depth 142 zkb_v = 2.0_wp * zk_v(ji,jj) * zbot_v ! 2k * bottom depth 133 143 ! 144 zke3_u = MAX(1.e-8_wp, 2.0_wp * zk_u(ji,jj) * e3u_n(ji,jj,jk)) ! 2k * thickness 145 zke3_v = MAX(1.e-8_wp, 2.0_wp * zk_v(ji,jj) * e3v_n(ji,jj,jk)) ! 2k * thickness 146 147 ! Depth attenuation .... do u component first.. 148 zdepth=zkb_u 149 zsqrt_depth = SQRT(zdepth) 150 zexp_depth = EXP(-zdepth) 151 zstokes_psi_u_bot = 1.0_wp - zexp_depth & 152 & - z_two_thirds*(zsqrtpi*zsqrt_depth*zdepth*ERFC(zsqrt_depth) & 153 & + 1.0_wp - (1.0_wp + zdepth)*zexp_depth) 154 zda_u = (zstokes_psi_u_bot - zstokes_psi_u_top(ji,jj))/zke3_u 155 zstokes_psi_u_top(ji,jj) = zstokes_psi_u_bot 156 157 ! ... and then v component 158 zdepth=zkb_v 159 zsqrt_depth = SQRT(zdepth) 160 zexp_depth = EXP(-zdepth) 161 zstokes_psi_v_bot = 1.0_wp - zexp_depth & 162 & - z_two_thirds*(zsqrtpi*zsqrt_depth*zdepth*ERFC(zsqrt_depth) & 163 & + 1.0_wp - (1.0_wp + zdepth)*zexp_depth) 164 zda_v = (zstokes_psi_v_bot - zstokes_psi_v_top(ji,jj))/zke3_v 165 zstokes_psi_v_top(ji,jj) = zstokes_psi_v_bot 166 ! 134 167 usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk) 135 168 vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk) … … 137 170 END DO 138 171 END DO 139 CALL lbc_lnk( usd(:,:,:), 'U', vsd(:,:,:), 'V', -1. ) 172 CALL lbc_lnk( usd(:,:,:), 'U', -1. ) 173 CALL lbc_lnk( vsd(:,:,:), 'V', -1. ) 174 175 140 176 ! 141 177 ! !== vertical Stokes Drift 3D velocity ==! … … 152 188 END DO 153 189 ! 154 IF( .NOT. AGRIF_Root() ) THEN 155 IF( nbondi == 1 .OR. nbondi == 2 ) ze3divh(nlci-1, : ,:) = 0._wp ! east 156 IF( nbondi == -1 .OR. nbondi == 2 ) ze3divh( 2 , : ,:) = 0._wp ! west 157 IF( nbondj == 1 .OR. nbondj == 2 ) ze3divh( : ,nlcj-1,:) = 0._wp ! north 158 IF( nbondj == -1 .OR. nbondj == 2 ) ze3divh( : , 2 ,:) = 0._wp ! south 159 ENDIF 190 #if defined key_agrif 191 IF( .NOT. Agrif_Root() ) THEN 192 IF( nbondi == -1 .OR. nbondi == 2 ) ze3divh( 2:nbghostcells+1,: ,:) = 0._wp ! west 193 IF( nbondi == 1 .OR. nbondi == 2 ) ze3divh( nlci-nbghostcells:nlci-1,:,:) = 0._wp ! east 194 IF( nbondj == -1 .OR. nbondj == 2 ) ze3divh( :,2:nbghostcells+1 ,:) = 0._wp ! south 195 IF( nbondj == 1 .OR. nbondj == 2 ) ze3divh( :,nlcj-nbghostcells:nlcj-1,:) = 0._wp ! north 196 ENDIF 197 #endif 160 198 ! 161 199 CALL lbc_lnk( ze3divh, 'T', 1. ) … … 185 223 CALL wrk_dealloc( jpi,jpj,jpk, ze3divh ) 186 224 CALL wrk_dealloc( jpi,jpj, zk_t, zk_u, zk_v, zu0_sd, zv0_sd ) 225 CALL wrk_dealloc( jpi,jpj, zstokes_psi_u_top, zstokes_psi_v_top) 187 226 ! 188 227 END SUBROUTINE sbc_stokes … … 227 266 ! 228 267 ! Read also wave number if needed, so that it is available in coupling routines 229 IF( ln_zdf qiao.AND. .NOT.cpl_wnum ) THEN268 IF( ln_zdfswm .AND. .NOT.cpl_wnum ) THEN 230 269 CALL fld_read( kt, nn_fsbc, sf_wn ) ! read wave parameters from external forcing 231 270 wnum(:,:) = sf_wn(1)%fnow(:,:,1) … … 342 381 ALLOCATE( div_sd(jpi,jpj) ) 343 382 ALLOCATE( tsd2d (jpi,jpj) ) 383 384 ut0sd(:,:) = 0._wp 385 vt0sd(:,:) = 0._wp 386 hsw(:,:) = 0._wp 387 wmp(:,:) = 0._wp 388 344 389 usd(:,:,:) = 0._wp 345 390 vsd(:,:,:) = 0._wp 346 391 wsd(:,:,:) = 0._wp 347 ! Wave number needed only if ln_zdf qiao=T392 ! Wave number needed only if ln_zdfswm=T 348 393 IF( .NOT. cpl_wnum ) THEN 349 394 ALLOCATE( sf_wn(1), STAT=ierror ) !* allocate and fill sf_wave with sn_wnum
Note: See TracChangeset
for help on using the changeset viewer.
