Changeset 15084 for NEMO/trunk/src/OCE/ISF/isfcavgam.F90
- Timestamp:
- 2021-07-05T21:48:42+02:00 (3 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/OCE/ISF/isfcavgam.F90 (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/ISF/isfcavgam.F90
r15082 r15084 14 14 USE isftbl , ONLY: isf_tbl 15 15 16 USE oce , ONLY: uu, vv , rn2 ! ocean dynamics and tracers16 USE oce , ONLY: uu, vv ! ocean dynamics 17 17 USE phycst , ONLY: grav, vkarmn ! physical constant 18 18 USE eosbn2 , ONLY: eos_rab ! equation of state … … 44 44 !!----------------------------------------------------------------------------------------------------- 45 45 ! 46 SUBROUTINE isfcav_gammats( Kmm, pttbl, pstbl, pqoce, pqfwf, p gt, pgs )46 SUBROUTINE isfcav_gammats( Kmm, pttbl, pstbl, pqoce, pqfwf, pRc, pgt, pgs ) 47 47 !!---------------------------------------------------------------------- 48 48 !! ** Purpose : compute the coefficient echange for heat and fwf flux … … 57 57 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pqoce, pqfwf ! isf heat and fwf 58 58 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pttbl, pstbl ! top boundary layer tracer 59 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pRc ! Richardson number 59 60 !!--------------------------------------------------------------------- 60 61 REAL(wp), DIMENSION(jpi,jpj) :: zutbl, zvtbl ! top boundary layer velocity … … 94 95 pgs(:,:) = rn_gammas0 95 96 CASE ( 'vel' ) ! gamma is proportional to u* 96 CALL gammats_vel ( zutbl, zvtbl, rCd0_top, r_ke0_top, pgt, pgs )97 CALL gammats_vel ( zutbl, zvtbl, rCd0_top, r_ke0_top, pgt, pgs ) 97 98 CASE ( 'vel_stab' ) ! gamma depends of stability of boundary layer and u* 98 CALL gammats_vel_stab (Kmm, pttbl, pstbl, zutbl, zvtbl, rCd0_top, r_ke0_top, pqoce, pqfwf, p gt, pgs )99 CALL gammats_vel_stab (Kmm, pttbl, pstbl, zutbl, zvtbl, rCd0_top, r_ke0_top, pqoce, pqfwf, pRc, pgt, pgs ) 99 100 CASE DEFAULT 100 101 CALL ctl_stop('STOP','method to compute gamma (cn_gammablk) is unknown (should not see this)') … … 153 154 END SUBROUTINE gammats_vel 154 155 155 SUBROUTINE gammats_vel_stab( Kmm, pttbl, pstbl, putbl, pvtbl, pCd, pke2, pqoce, pqfwf, & ! <<== in156 & pgt , pgs ) ! ==>> out gammats [m/s]156 SUBROUTINE gammats_vel_stab( Kmm, pttbl, pstbl, putbl, pvtbl, pCd, pke2, pqoce, pqfwf, pRc, & ! <<== in 157 & pgt , pgs ) ! ==>> out gammats [m/s] 157 158 !!---------------------------------------------------------------------- 158 159 !! ** Purpose : compute the coefficient echange coefficient … … 171 172 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: putbl, pvtbl ! velocity in the losch top boundary layer 172 173 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pttbl, pstbl ! tracer in the losch top boundary layer 174 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pRc ! Richardson number 173 175 !!--------------------------------------------------------------------- 174 176 INTEGER :: ji, jj ! loop index 175 177 INTEGER :: ikt ! local integer 176 178 REAL(wp) :: zdku, zdkv ! U, V shear 177 REAL(wp) :: zPr, zSc , zRc ! Prandtl, Scmidth and Richardsonnumber179 REAL(wp) :: zPr, zSc ! Prandtl and Scmidth number 178 180 REAL(wp) :: zmob, zmols ! Monin Obukov length, coriolis factor at T point 179 181 REAL(wp) :: zbuofdep, zhnu ! Bouyancy length scale, sublayer tickness … … 190 192 !!--------------------------------------------------------------------- 191 193 ! 192 ! compute ustar193 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )194 zustar(ji,jj) = SQRT( pCd(ji,jj) * ( putbl(ji,jj) * putbl(ji,jj) + pvtbl(ji,jj) * pvtbl(ji,jj) + pke2 ) )195 END_2D196 !197 ! output ustar198 CALL iom_put('isfustar',zustar(:,:))199 !200 194 ! compute Pr and Sc number (eq ??) 201 195 zPr = 13.8_wp … … 207 201 ! 208 202 ! compute gamma 209 DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls ) 203 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 204 210 205 ikt = mikt(ji,jj) 206 207 ! compute ustar 208 zustar(ji,jj) = SQRT( pCd(ji,jj) * ( putbl(ji,jj) * putbl(ji,jj) + pvtbl(ji,jj) * pvtbl(ji,jj) + pke2 ) ) 211 209 212 210 IF( zustar(ji,jj) == 0._wp ) THEN ! only for kt = 1 I think … … 214 212 pgs(ji,jj) = rn_gammas0 215 213 ELSE 216 ! compute Rc number (as done in zdfric.F90)217 !!gm better to do it like in the new zdfric.F90 i.e. avm weighted Ri computation218 zcoef = 0.5_wp / e3w(ji,jj,ikt+1,Kmm)219 ! ! shear of horizontal velocity220 zdku = zcoef * ( uu(ji-1,jj ,ikt ,Kmm) + uu(ji,jj,ikt ,Kmm) &221 & -uu(ji-1,jj ,ikt+1,Kmm) - uu(ji,jj,ikt+1,Kmm) )222 zdkv = zcoef * ( vv(ji ,jj-1,ikt ,Kmm) + vv(ji,jj,ikt ,Kmm) &223 & -vv(ji ,jj-1,ikt+1,Kmm) - vv(ji,jj,ikt+1,Kmm) )224 ! ! richardson number (minimum value set to zero)225 zRc = MAX(rn2(ji,jj,ikt+1), 0._wp) / MAX( zdku*zdku + zdkv*zdkv, zeps )226 227 214 ! compute bouyancy 228 215 zts(jp_tem) = pttbl(ji,jj) … … 244 231 ! 245 232 ! compute eta* (stability parameter) (Eq ??) 246 zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff_f(ji,jj)) * zmols * zRc ), 0._wp))) 233 zetastar = 1._wp / ( SQRT(1._wp + MAX( 0._wp, zxsiN * zustar(ji,jj) & 234 & / MAX( 1.e-20, ABS(ff_t(ji,jj)) * zmols * pRc(ji,jj) ) ))) 247 235 ! 248 236 ! compute the sublayer thickness (Eq ??) 249 zhnu = 5 * znu / zustar(ji,jj)237 zhnu = 5 * znu / MAX( 1.e-20, zustar(ji,jj) ) 250 238 ! 251 239 ! compute gamma turb (Eq ??) 252 zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff_f(ji,jj)) * zhnu )) &240 zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / MAX( 1.e-10, ABS(ff_t(ji,jj)) * zhnu )) & 253 241 & + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn 254 242 ! … … 258 246 END IF 259 247 END_2D 248 ! output ustar 249 CALL iom_put('isfustar',zustar(:,:)) 260 250 261 251 END SUBROUTINE gammats_vel_stab
Note: See TracChangeset
for help on using the changeset viewer.
