Changeset 15082
- Timestamp:
- 2021-07-05T18:25:26+02:00 (3 years ago)
- Location:
- NEMO/trunk/src/OCE/ISF
- Files:
-
- 2 edited
-
isfcav.F90 (modified) (4 diffs)
-
isfcavgam.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/ISF/isfcav.F90
r15004 r15082 74 74 REAL(wp) :: zerr 75 75 REAL(wp), DIMENSION(jpi,jpj) :: zqlat, zqoce, zqhc, zqh ! heat fluxes 76 REAL(wp), DIMENSION(jpi,jpj) :: zq oce_b!76 REAL(wp), DIMENSION(jpi,jpj) :: zqh_b ! 77 77 REAL(wp), DIMENSION(jpi,jpj) :: zgammat, zgammas ! exchange coeficient 78 78 REAL(wp), DIMENSION(jpi,jpj) :: zttbl, zstbl ! temp. and sal. in top boundary layer … … 88 88 ! 89 89 ! initialisation 90 IF (TRIM(cn_gammablk) == 'vel_stab' ) zqoce_b (:,:) = ptsc(:,:,jp_tem) * rho0_rcp ! last time step total heat fluxes (to speed up convergence) 90 IF ( TRIM(cn_gammablk) == 'vel_stab' ) THEN 91 zqoce(:,:) = -pqfwf(:,:) * rLfusisf ! 92 zqh_b(:,:) = ptsc(:,:,jp_tem) * rho0_rcp ! last time step total heat fluxes (to speed up convergence) 93 ENDIF 91 94 ! 92 95 ! compute ice shelf melting … … 112 115 CASE ( 'vel_stab' ) 113 116 ! compute error between 2 iterations 114 zerr = MAXVAL(ABS(zq oce(:,:) - zqoce_b(:,:)))117 zerr = MAXVAL(ABS(zqhc(:,:)+zqoce(:,:) - zqh_b(:,:))) 115 118 ! 116 119 ! define if iteration needed … … 121 124 ELSE ! converge is not yet achieve 122 125 nit = nit + 1 123 zq oce_b(:,:) =zqoce(:,:)126 zqh_b(:,:) = zqhc(:,:)+zqoce(:,:) 124 127 END IF 125 128 END SELECT -
NEMO/trunk/src/OCE/ISF/isfcavgam.F90
r13237 r15082 30 30 PUBLIC isfcav_gammats 31 31 32 !! * Substitutions 33 # include "do_loop_substitute.h90" 32 34 # include "domzgr_substitute.h90" 33 35 !!---------------------------------------------------------------------- … … 133 135 REAL(wp), INTENT(in ) :: pke2 ! background velocity 134 136 !!--------------------------------------------------------------------- 137 INTEGER :: ji, jj ! loop index 135 138 REAL(wp), DIMENSION(jpi,jpj) :: zustar 136 139 !!--------------------------------------------------------------------- 137 140 ! 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 141 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 142 ! compute ustar (AD15 eq. 27) 143 zustar(ji,jj) = SQRT( pCd(ji,jj) * ( putbl(ji,jj) * putbl(ji,jj) + pvtbl(ji,jj) * pvtbl(ji,jj) + pke2 ) ) * mskisf_cav(ji,jj) 144 ! 145 ! Compute gammats 146 pgt(ji,jj) = zustar(ji,jj) * rn_gammat0 147 pgs(ji,jj) = zustar(ji,jj) * rn_gammas0 148 END_2D 144 149 ! 145 150 ! output ustar … … 186 191 ! 187 192 ! compute ustar 188 zustar(:,:) = SQRT( pCd * ( putbl(:,:) * putbl(:,:) + pvtbl(:,:) * pvtbl(:,:) + pke2 ) ) 193 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_2D 189 196 ! 190 197 ! output ustar … … 200 207 ! 201 208 ! 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) 209 DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls ) 210 ikt = mikt(ji,jj) 211 212 IF( zustar(ji,jj) == 0._wp ) THEN ! only for kt = 1 I think 213 pgt(ji,jj) = rn_gammat0 214 pgs(ji,jj) = rn_gammas0 215 ELSE 216 ! compute Rc number (as done in zdfric.F90) 211 217 !!gm better to do it like in the new zdfric.F90 i.e. avm weighted Ri computation 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 218 zcoef = 0.5_wp / e3w(ji,jj,ikt+1,Kmm) 219 ! ! shear of horizontal velocity 220 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 ! compute bouyancy 228 zts(jp_tem) = pttbl(ji,jj) 229 zts(jp_sal) = pstbl(ji,jj) 230 zdep = gdepw(ji,jj,ikt,Kmm) 231 ! 232 CALL eos_rab( zts, zdep, zab, Kmm ) 233 ! 234 ! compute length scale (Eq ??) 235 zbuofdep = grav * ( zab(jp_tem) * pqoce(ji,jj) - zab(jp_sal) * pqfwf(ji,jj) ) 236 ! 237 ! compute Monin Obukov Length 238 ! Maximum boundary layer depth (Eq ??) 239 zhmax = gdept(ji,jj,mbkt(ji,jj),Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm) -0.001_wp 240 ! 241 ! Compute Monin obukhov length scale at the surface and Ekman depth: (Eq ??) 242 zmob = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps)) 243 zmols = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt) 244 ! 245 ! 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))) 247 ! 248 ! compute the sublayer thickness (Eq ??) 249 zhnu = 5 * znu / zustar(ji,jj) 250 ! 251 ! compute gamma turb (Eq ??) 252 zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff_f(ji,jj)) * zhnu )) & 247 253 & + 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 254 ! 255 ! compute gammats 256 pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet) 257 pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles) 258 END IF 259 END_2D 255 260 256 261 END SUBROUTINE gammats_vel_stab
Note: See TracChangeset
for help on using the changeset viewer.
