Changeset 15109
- Timestamp:
- 2021-07-08T14:56:28+02:00 (3 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/tests/ISOMIP+/MY_SRC/isfcavgam.F90 (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/tests/ISOMIP+/MY_SRC/isfcavgam.F90
r13583 r15109 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 … … 30 30 PUBLIC isfcav_gammats 31 31 32 !! * Substitutions 33 # include "do_loop_substitute.h90" 32 34 # include "domzgr_substitute.h90" 33 35 !!---------------------------------------------------------------------- … … 42 44 !!----------------------------------------------------------------------------------------------------- 43 45 ! 44 SUBROUTINE isfcav_gammats( Kmm, pttbl, pstbl, pqoce, pqfwf, p gt, pgs )46 SUBROUTINE isfcav_gammats( Kmm, pttbl, pstbl, pqoce, pqfwf, pRc, pgt, pgs ) 45 47 !!---------------------------------------------------------------------- 46 48 !! ** Purpose : compute the coefficient echange for heat and fwf flux … … 55 57 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pqoce, pqfwf ! isf heat and fwf 56 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 57 60 !!--------------------------------------------------------------------- 58 61 REAL(wp), DIMENSION(jpi,jpj) :: zutbl, zvtbl ! top boundary layer velocity … … 92 95 pgs(:,:) = rn_gammas0 93 96 CASE ( 'vel' ) ! gamma is proportional to u* 94 CALL gammats_vel ( zutbl, zvtbl, rCd0_top, rn_vtide**2, pgt, pgs )97 CALL gammats_vel ( zutbl, zvtbl, rCd0_top, rn_vtide**2, pgt, pgs ) 95 98 CASE ( 'vel_stab' ) ! gamma depends of stability of boundary layer and u* 96 CALL gammats_vel_stab (Kmm, pttbl, pstbl, zutbl, zvtbl, rCd0_top, rn_vtide**2, pqoce, pqfwf, p gt, pgs )99 CALL gammats_vel_stab (Kmm, pttbl, pstbl, zutbl, zvtbl, rCd0_top, rn_vtide**2, pqoce, pqfwf, pRc, pgt, pgs ) 97 100 CASE DEFAULT 98 101 CALL ctl_stop('STOP','method to compute gamma (cn_gammablk) is unknown (should not see this)') … … 133 136 REAL(wp), INTENT(in ) :: pke2 ! background velocity 134 137 !!--------------------------------------------------------------------- 138 INTEGER :: ji, jj ! loop index 135 139 REAL(wp), DIMENSION(jpi,jpj) :: zustar 136 140 !!--------------------------------------------------------------------- 137 141 ! 138 ! compute ustar (AD15 eq. 27) 139 zustar(:,:) = SQRT( pCd(:,:) * ( putbl(:,:) * putbl(:,:) + pvtbl(:,:) * pvtbl(:,:) + pke2 ) ) * mskisf_cav(:,:) 140 ! 141 ! Compute gammats 142 pgt(:,:) = zustar(:,:) * rn_gammat0 143 pgs(:,:) = zustar(:,:) * rn_gammas0 142 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 143 ! compute ustar (AD15 eq. 27) 144 zustar(ji,jj) = SQRT( pCd(ji,jj) * ( putbl(ji,jj) * putbl(ji,jj) + pvtbl(ji,jj) * pvtbl(ji,jj) + pke2 ) ) * mskisf_cav(ji,jj) 145 ! 146 ! Compute gammats 147 pgt(ji,jj) = zustar(ji,jj) * rn_gammat0 148 pgs(ji,jj) = zustar(ji,jj) * rn_gammas0 149 END_2D 144 150 ! 145 151 ! output ustar … … 148 154 END SUBROUTINE gammats_vel 149 155 150 SUBROUTINE gammats_vel_stab( Kmm, pttbl, pstbl, putbl, pvtbl, pCd, pke2, pqoce, pqfwf, & ! <<== in151 & 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] 152 158 !!---------------------------------------------------------------------- 153 159 !! ** Purpose : compute the coefficient echange coefficient … … 166 172 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: putbl, pvtbl ! velocity in the losch top boundary layer 167 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 168 175 !!--------------------------------------------------------------------- 169 176 INTEGER :: ji, jj ! loop index 170 177 INTEGER :: ikt ! local integer 171 178 REAL(wp) :: zdku, zdkv ! U, V shear 172 REAL(wp) :: zPr, zSc , zRc ! Prandtl, Scmidth and Richardsonnumber179 REAL(wp) :: zPr, zSc ! Prandtl and Scmidth number 173 180 REAL(wp) :: zmob, zmols ! Monin Obukov length, coriolis factor at T point 174 181 REAL(wp) :: zbuofdep, zhnu ! Bouyancy length scale, sublayer tickness … … 185 192 !!--------------------------------------------------------------------- 186 193 ! 187 ! compute ustar188 zustar(:,:) = SQRT( pCd * ( putbl(:,:) * putbl(:,:) + pvtbl(:,:) * pvtbl(:,:) + pke2 ) )189 !190 ! output ustar191 CALL iom_put('isfustar',zustar(:,:))192 !193 194 ! compute Pr and Sc number (eq ??) 194 195 zPr = 13.8_wp … … 200 201 ! 201 202 ! compute gamma 202 DO ji = 2, jpi 203 DO jj = 2, jpj 204 ikt = mikt(ji,jj) 205 206 IF( zustar(ji,jj) == 0._wp ) THEN ! only for kt = 1 I think 207 pgt = rn_gammat0 208 pgs = rn_gammas0 209 ELSE 210 ! compute Rc number (as done in zdfric.F90) 211 !!gm better to do it like in the new zdfric.F90 i.e. avm weighted Ri computation 212 zcoef = 0.5_wp / e3w(ji,jj,ikt+1,Kmm) 213 ! ! shear of horizontal velocity 214 zdku = zcoef * ( uu(ji-1,jj ,ikt ,Kmm) + uu(ji,jj,ikt ,Kmm) & 215 & -uu(ji-1,jj ,ikt+1,Kmm) - uu(ji,jj,ikt+1,Kmm) ) 216 zdkv = zcoef * ( vv(ji ,jj-1,ikt ,Kmm) + vv(ji,jj,ikt ,Kmm) & 217 & -vv(ji ,jj-1,ikt+1,Kmm) - vv(ji,jj,ikt+1,Kmm) ) 218 ! ! richardson number (minimum value set to zero) 219 zRc = MAX(rn2(ji,jj,ikt+1), 0._wp) / MAX( zdku*zdku + zdkv*zdkv, zeps ) 220 221 ! compute bouyancy 222 zts(jp_tem) = pttbl(ji,jj) 223 zts(jp_sal) = pstbl(ji,jj) 224 zdep = gdepw(ji,jj,ikt,Kmm) 225 ! 226 CALL eos_rab( zts, zdep, zab, Kmm ) 227 ! 228 ! compute length scale (Eq ??) 229 zbuofdep = grav * ( zab(jp_tem) * pqoce(ji,jj) - zab(jp_sal) * pqfwf(ji,jj) ) 230 ! 231 ! compute Monin Obukov Length 232 ! Maximum boundary layer depth (Eq ??) 233 zhmax = gdept(ji,jj,mbkt(ji,jj),Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm) -0.001_wp 234 ! 235 ! Compute Monin obukhov length scale at the surface and Ekman depth: (Eq ??) 236 zmob = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps)) 237 zmols = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt) 238 ! 239 ! compute eta* (stability parameter) (Eq ??) 240 zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff_f(ji,jj)) * zmols * zRc ), 0._wp))) 241 ! 242 ! compute the sublayer thickness (Eq ??) 243 zhnu = 5 * znu / zustar(ji,jj) 244 ! 245 ! compute gamma turb (Eq ??) 246 zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff_f(ji,jj)) * zhnu )) & 203 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 204 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 ) ) 209 210 IF( zustar(ji,jj) == 0._wp ) THEN ! only for kt = 1 I think 211 pgt(ji,jj) = rn_gammat0 212 pgs(ji,jj) = rn_gammas0 213 ELSE 214 ! compute bouyancy 215 zts(jp_tem) = pttbl(ji,jj) 216 zts(jp_sal) = pstbl(ji,jj) 217 zdep = gdepw(ji,jj,ikt,Kmm) 218 ! 219 CALL eos_rab( zts, zdep, zab, Kmm ) 220 ! 221 ! compute length scale (Eq ??) 222 zbuofdep = grav * ( zab(jp_tem) * pqoce(ji,jj) - zab(jp_sal) * pqfwf(ji,jj) ) 223 ! 224 ! compute Monin Obukov Length 225 ! Maximum boundary layer depth (Eq ??) 226 zhmax = gdept(ji,jj,mbkt(ji,jj),Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm) -0.001_wp 227 ! 228 ! Compute Monin obukhov length scale at the surface and Ekman depth: (Eq ??) 229 zmob = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps)) 230 zmols = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt) 231 ! 232 ! compute eta* (stability parameter) (Eq ??) 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) ) ))) 235 ! 236 ! compute the sublayer thickness (Eq ??) 237 zhnu = 5 * znu / MAX( 1.e-20, zustar(ji,jj) ) 238 ! 239 ! compute gamma turb (Eq ??) 240 zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / MAX( 1.e-10, ABS(ff_t(ji,jj)) * zhnu )) & 247 241 & + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn 248 ! 249 ! compute gammats 250 pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet) 251 pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles) 252 END IF 253 END DO 254 END DO 242 ! 243 ! compute gammats 244 pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet) 245 pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles) 246 END IF 247 END_2D 248 ! output ustar 249 CALL iom_put('isfustar',zustar(:,:)) 255 250 256 251 END SUBROUTINE gammats_vel_stab
Note: See TracChangeset
for help on using the changeset viewer.
