- Timestamp:
- 2016-10-25T18:15:30+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/r5518_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
r5215 r7099 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 8 !!---------------------------------------------------------------------- 9 USE iom ! I/O manager library 10 USE in_out_manager ! I/O manager 11 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 18 !!---------------------------------------------------------------------- 19 !! sbc_wave : read drag coefficient from wave model in netcdf files 20 !!---------------------------------------------------------------------- 6 !! History : 3.3 ! 2011-09 (Adani M) Original code: Drag Coefficient 7 !! : 3.4 ! 2012-10 (Adani M) Stokes Drift 8 !! 3.6 ! 2014-09 (Clementi E, Oddo P)New Stokes Drift Computation 9 !!---------------------------------------------------------------------- 10 11 !!---------------------------------------------------------------------- 12 !! sbc_wave : wave data from wave model in netcdf files 13 !!---------------------------------------------------------------------- 14 USE oce ! 15 USE sbc_oce ! Surface boundary condition: ocean fields 16 USE bdy_oce ! 17 USE domvvl ! 18 ! 19 USE iom ! I/O manager library 20 USE in_out_manager ! I/O manager 21 USE lib_mpp ! distribued memory computing library 22 USE fldread ! read input fields 23 USE wrk_nemo ! 24 USE phycst ! physical constants 21 25 22 26 IMPLICIT NONE 23 27 PRIVATE 24 28 25 PUBLIC sbc_wave ! routine called in sbc _blk_core or sbc_blk_mfs29 PUBLIC sbc_wave ! routine called in sbcmod 26 30 27 INTEGER , PARAMETER :: jpfld = 3 ! maximum number of files to read for srokes drift31 INTEGER , PARAMETER :: jpfld = 4 ! number of files to read for stokes drift 28 32 INTEGER , PARAMETER :: jp_usd = 1 ! index of stokes drift (i-component) (m/s) at T-point 29 33 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 34 INTEGER , PARAMETER :: jp_swh = 3 ! index of significant wave hight (m) at T-point 35 INTEGER , PARAMETER :: jp_wmp = 4 ! index of mean wave period (s) at T-point 36 37 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_cd ! structure of input fields (file informations, fields read) Drag Coefficient 38 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_sd ! structure of input fields (file informations, fields read) Stokes Drift 39 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_wn ! structure of input fields (file informations, fields read) wave number for Qiao 40 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_tauoc ! structure of input fields (file informations, fields read) normalized wave stress into the ocean 41 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: cdn_wave 42 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: swh,wmp, wnum 43 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: tauoc_wave 44 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: tsd2d 45 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: usd2d, vsd2d 46 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: usd3d, vsd3d, wsd3d 47 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: zusd2dt, zvsd2dt 36 48 37 49 !! * Substitutions 38 50 # include "domzgr_substitute.h90" 39 !!---------------------------------------------------------------------- 40 !! NEMO/OPA 4.0 , NEMO Consortium (2011) 51 # include "vectopt_loop_substitute.h90" 52 !!---------------------------------------------------------------------- 53 !! NEMO/OPA 3.7 , NEMO Consortium (2014) 41 54 !! $Id$ 42 55 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 46 59 SUBROUTINE sbc_wave( kt ) 47 60 !!--------------------------------------------------------------------- 48 !! *** ROUTINE sbc_ apr***61 !! *** ROUTINE sbc_wave *** 49 62 !! 50 !! ** Purpose : read drag coefficientfrom wave model in netcdf files.63 !! ** Purpose : read wave parameters from wave model in netcdf files. 51 64 !! 52 65 !! ** Method : - Read namelist namsbc_wave 53 66 !! - Read Cd_n10 fields in netcdf files 54 67 !! - Read stokes drift 2d in netcdf files 55 !! - Read wave number 56 !! - Compute 3d stokes drift using monochromatic57 !! ** action :58 !! 68 !! - Read wave number in netcdf files 69 !! - Compute 3d stokes drift using Breivik et al.,2014 70 !! formulation 71 !! ** action 59 72 !!--------------------------------------------------------------------- 60 USE oce, ONLY : un,vn,hdivn,rotn61 USE divcur 62 USE wrk_nemo63 #if defined key_bdy 64 USE bdy_oce, ONLY : bdytmask65 #endif 66 INTEGER , INTENT( in ) :: kt ! ocean time step67 INTEGER :: ierror ! return error code68 INTEGER :: ifpr, jj,ji,jk69 INTEGER :: ios ! Local integer output status for namelist read70 REAL(wp), DIMENSION(:,:,:),POINTER :: udummy,vdummy,hdivdummy,rotdummy71 REAL :: z2dt,z1_2dt73 USE zdf_oce, ONLY : ln_zdfqiao 74 75 INTEGER, INTENT( in ) :: kt ! ocean time step 76 ! 77 INTEGER :: ierror ! return error code 78 INTEGER :: ifpr, jj,ji,jk 79 INTEGER :: ios ! Local integer output status for namelist read 80 ! 81 CHARACTER(len=100) :: cn_dir ! Root directory for location of drag coefficient files 82 REAL(wp) :: ztransp, zsp0, zk, zus, zvs 83 REAL(wp), DIMENSION(jpi,jpj) :: zfac 84 REAL(wp), DIMENSION(:,:,:), POINTER :: ze3hdiv ! 3D workspace 72 85 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 86 TYPE(FLD_N) :: sn_cdg, sn_usd, sn_vsd, & 87 & sn_swh, sn_wmp, sn_wnum, sn_tauoc ! informations about the fields to be read 88 !! 89 NAMELIST/namsbc_wave/ sn_cdg, cn_dir, sn_usd, sn_vsd, sn_swh, sn_wmp, sn_wnum, sn_tauoc 75 90 !!--------------------------------------------------------------------- 76 NAMELIST/namsbc_wave/ sn_cdg, cn_dir, sn_usd, sn_vsd, sn_wn77 !!---------------------------------------------------------------------78 79 !!----------------------------------------------------------------------80 !81 91 ! 82 92 ! ! -------------------- ! … … 92 102 IF(lwm) WRITE ( numond, namsbc_wave ) 93 103 ! 94 95 104 IF ( ln_cdgw ) THEN 96 105 ALLOCATE( sf_cd(1), STAT=ierror ) !* allocate and fill sf_wave with sn_cdg … … 102 111 ALLOCATE( cdn_wave(jpi,jpj) ) 103 112 cdn_wave(:,:) = 0.0 113 ENDIF 114 115 IF ( ln_tauoc ) THEN 116 ALLOCATE( sf_tauoc(1), STAT=ierror ) !* allocate and fill sf_wave with sn_tauoc 117 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 118 ! 119 ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1) ) 120 IF( sn_tauoc%ln_tint ) ALLOCATE( sf_tauoc(1)%fdta(jpi,jpj,1,2) ) 121 CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' ) 122 ALLOCATE( tauoc_wave(jpi,jpj) ) 123 tauoc_wave(:,:) = 0.0 104 124 ENDIF 125 105 126 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 127 slf_i(jp_usd) = sn_usd ; slf_i(jp_vsd) = sn_vsd; 128 slf_i(jp_swh) = sn_swh ; slf_i(jp_wmp) = sn_wmp; 129 ALLOCATE( sf_sd(jpfld), STAT=ierror ) !* allocate and fill sf_sd with stokes drift 108 130 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 109 131 ! … … 112 134 IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) ) 113 135 END DO 136 114 137 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))138 ALLOCATE( usd2d(jpi,jpj),vsd2d(jpi,jpj) ) 116 139 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 119 ENDIF 120 ENDIF 121 ! 122 ! 123 IF ( ln_cdgw ) THEN 124 CALL fld_read( kt, nn_fsbc, sf_cd ) !* read drag coefficient from external forcing 140 ALLOCATE( swh(jpi,jpj), wmp(jpi,jpj) ) 141 ALLOCATE( zusd2dt(jpi,jpj), zvsd2dt(jpi,jpj) ) 142 usd3d(:,:,:) = 0._wp ; usd2d(:,:) = 0._wp ; 143 vsd3d(:,:,:) = 0._wp ; vsd2d(:,:) = 0._wp ; 144 wsd3d(:,:,:) = 0._wp ; 145 swh (:,:) = 0._wp ; wmp (:,:) = 0._wp ; 146 IF ( ln_zdfqiao ) THEN !== Vertical mixing enhancement using Qiao,2010 ==! 147 ALLOCATE( sf_wn(1), STAT=ierror ) !* allocate and fill sf_wave with sn_wnum 148 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave: unable toallocate sf_wave structure' ) 149 ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1) ) 150 IF( sn_wnum%ln_tint ) ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) ) 151 CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' ) 152 ALLOCATE( wnum(jpi,jpj),tsd2d(jpi,jpj) ) 153 wnum(:,:) = 0._wp ; tsd2d(:,:) = 0._wp 154 ENDIF 155 ENDIF 156 ENDIF 157 ! 158 IF ( ln_cdgw ) THEN !== Neutral drag coefficient ==! 159 CALL fld_read( kt, nn_fsbc, sf_cd ) ! read from external forcing 125 160 cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1) 126 161 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) ) 162 163 IF ( ln_tauoc ) THEN !== Wave induced stress ==! 164 CALL fld_read( kt, nn_fsbc, sf_tauoc ) !* read wave norm stress from external forcing 165 tauoc_wave(:,:) = sf_tauoc(1)%fnow(:,:,1) 166 ENDIF 167 168 IF ( ln_sdw ) THEN !== Computation of the 3d Stokes Drift ==! 169 ! 170 CALL fld_read( kt, nn_fsbc, sf_sd ) !* read wave parameters from external forcing 171 swh(:,:) = sf_sd(jp_swh)%fnow(:,:,1) ! significant wave height 172 wmp(:,:) = sf_sd(jp_wmp)%fnow(:,:,1) ! wave mean period 173 zusd2dt(:,:) = sf_sd(jp_usd)%fnow(:,:,1) ! 2D zonal Stokes Drift at T point 174 zvsd2dt(:,:) = sf_sd(jp_vsd)%fnow(:,:,1) ! 2D meridional Stokes Drift at T point 175 ! 176 !== Computation of the 3d Stokes Drift according to Breivik et al.,2014 177 !(DOI: 10.1175/JPO-D-14-0020.1)==! 178 ! 179 CALL wrk_alloc( jpi,jpj,jpk, ze3hdiv ) 180 DO jk = 1, jpk 181 DO jj = 1, jpj 182 DO ji = 1, jpi 183 ! On T grid 184 ! Stokes transport speed estimated from Hs and Tmean 185 ztransp = 2.0_wp*rpi*swh(ji,jj)**2.0_wp/(16.0_wp*MAX(wmp(ji,jj),0.0000001_wp)) 186 ! Stokes surface speed 187 zsp0 = SQRT( sf_sd(jp_usd)%fnow(ji,jj,1)**2 + sf_sd(jp_vsd)%fnow(ji,jj,1)**2) 188 ! Wavenumber scale 189 zk = ABS(zsp0)/MAX(ABS(5.97_wp*ztransp),0.0000001_wp) 190 ! Depth attenuation 191 zfac(ji,jj) = EXP(-2.0_wp*zk*fsdept(ji,jj,jk))/(1.0_wp+8.0_wp*zk*fsdept(ji,jj,jk)) 192 END DO 193 END DO 194 195 DO jj = 1, jpjm1 196 DO ji = 1, jpim1 197 ! Into the U and V Grid 198 zus = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zfac(ji,jj) * tmask(ji,jj,1) & 199 & + zfac(ji+1,jj) * tmask(ji+1,jj,1) ) 200 201 zvs = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zfac(ji,jj) * tmask(ji,jj,1) & 202 & + zfac(ji,jj+1) * tmask(ji,jj+1,1) ) 203 204 usd2d(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zusd2dt(ji,jj) * tmask(ji,jj,1) & 205 & + zusd2dt(ji+1,jj) * tmask(ji+1,jj,1) ) 206 207 vsd2d(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zvsd2dt(ji,jj) * tmask(ji,jj,1) & 208 & + zvsd2dt(ji,jj+1) * tmask(ji,jj+1,1) ) 209 210 usd3d(ji,jj,jk) = usd2d(ji,jj) * zus 211 vsd3d(ji,jj,jk) = vsd2d(ji,jj) * zvs 212 END DO 146 213 END DO 147 214 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 215 ! 216 CALL lbc_lnk( usd3d(:,:,:), 'U', -1. ) 217 CALL lbc_lnk( vsd3d(:,:,:), 'V', -1. ) 218 ! 219 DO jk = 1, jpkm1 ! e3t * Horizontal divergence 220 DO jj = 2, jpjm1 221 DO ji = fs_2, fs_jpim1 ! vector opt. 222 ze3hdiv(ji,jj,jk) = ( e2u(ji ,jj) * fse3u_n(ji ,jj,jk) * usd3d(ji ,jj,jk) & 223 & - e2u(ji-1,jj) * fse3u_n(ji-1,jj,jk) * usd3d(ji-1,jj,jk) & 224 & + e1v(ji,jj ) * fse3v_n(ji,jj ,jk) * vsd3d(ji,jj ,jk) & 225 & - e1v(ji,jj-1) * fse3v_n(ji,jj-1,jk) * vsd3d(ji,jj-1,jk) ) * r1_e1e2t(ji,jj) 226 END DO 227 END DO 228 IF( .NOT. AGRIF_Root() ) THEN 229 IF( nbondi == 1 .OR. nbondi == 2 ) ze3hdiv(nlci-1, : ,jk) = 0._wp ! east 230 IF( nbondi == -1 .OR. nbondi == 2 ) ze3hdiv( 2 , : ,jk) = 0._wp ! west 231 IF( nbondj == 1 .OR. nbondj == 2 ) ze3hdiv( : ,nlcj-1,jk) = 0._wp ! north 232 IF( nbondj == -1 .OR. nbondj == 2 ) ze3hdiv( : , 2 ,jk) = 0._wp ! south 233 ENDIF 234 END DO 235 CALL lbc_lnk( ze3hdiv, 'T', 1. ) 236 ! 237 DO jk = jpkm1, 1, -1 !* integrate from the bottom the e3t * hor. divergence 238 wsd3d(:,:,jk) = wsd3d(:,:,jk+1) - ze3hdiv(:,:,jk) 239 END DO 181 240 #if defined key_bdy 182 wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:) 241 IF( lk_bdy ) THEN 242 DO jk = 1, jpkm1 243 wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:) 244 END DO 245 ENDIF 183 246 #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 247 CALL wrk_dealloc( jpi,jpj,jpk, ze3hdiv ) 248 249 IF ( ln_zdfqiao ) THEN 250 wnum(:,:) = sf_wn(1)%fnow(:,:,1) 251 252 ! Calculate the module of the stokes drift on T grid 253 !------------------------------------------------- 254 DO jj = 1, jpj 255 DO ji = 1, jpi 256 tsd2d(ji,jj) = ((sf_sd(jp_usd)%fnow(ji,jj,1) * tmask(ji,jj,1))**2.0 + & 257 & (sf_sd(jp_vsd)%fnow(ji,jj,1) * tmask(ji,jj,1))**2.0)**0.5 258 END DO 259 END DO 260 ENDIF 261 ! 262 ENDIF 263 ! 191 264 END SUBROUTINE sbc_wave 192 265
Note: See TracChangeset
for help on using the changeset viewer.
