Changeset 15515
- Timestamp:
- 2021-11-16T14:15:05+01:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_vp.F90
r15499 r15515 164 164 ! 165 165 REAL(wp), DIMENSION(jpi,jpj) :: zdeltat, zdeltastar_t ! Delta & Delta* at T-points 166 REAL(wp), DIMENSION(jpi,jpj) :: zten _i ! Tension166 REAL(wp), DIMENSION(jpi,jpj) :: ztens, zshear ! Tension, shear 167 167 REAL(wp), DIMENSION(jpi,jpj) :: zp_deltastar_t ! P/delta* at T points 168 168 REAL(wp), DIMENSION(jpi,jpj) :: zzt, zet ! Viscosity pre-factors at T points … … 538 538 539 539 ! Temporary zef factor at F-point 540 zef(ji,jj) = zp_deltastar_f * r1_e1e2f(ji,jj) * z1_ecc2 * zfmask(ji,jj) 540 zef(ji,jj) = zp_deltastar_f * r1_e1e2f(ji,jj) * z1_ecc2 * zfmask(ji,jj) * 0.5_wp ! BUG missing un facteur 2 Nov 12 541 541 542 542 END DO … … 1219 1219 zdt2 = zdt * zdt 1220 1220 1221 zten _i(ji,jj)= zdt1221 ztens(ji,jj) = zdt 1222 1222 1223 1223 ! shear**2 at T points (doc eq. A16) … … 1226 1226 & ) * 0.25_wp * r1_e1e2t(ji,jj) 1227 1227 1228 ! shear at T points 1229 pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) 1228 ! maximum shear rate at T points (includees tension, output only) 1229 pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) ! i think this is maximum shear rate and not actual shear --- i'm not totally sure here 1230 1231 ! shear at T-points 1232 zshear(ji,jj) = SQRT( zds2 ) 1230 1233 1231 1234 ! divergence at T points … … 1233 1236 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 1234 1237 & ) * r1_e1e2t(ji,jj) 1238 1239 zdiv2 = pdivu_i(ji,jj) * pdivu_i(ji,jj) 1235 1240 1236 1241 ! delta at T points 1237 zdelta = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )1242 zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 1238 1243 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0 1239 1244 … … 1261 1266 ! 1262 1267 ! ---- Sea ice stresses at T-points 1263 IF ( iom_use('normstr') .OR. iom_use('sheastr') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN 1268 IF ( iom_use('normstr') .OR. iom_use('sheastr') .OR. & 1269 & iom_use('intstrx') .OR. iom_use('intstry') .OR. & 1270 & iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 1264 1271 1265 1272 DO jj = 2, jpj - 1 … … 1268 1275 zfac = zp_deltastar_t(ji,jj) 1269 1276 zs1(ji,jj) = zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 1270 zs2(ji,jj) = zfac * z1_ecc2 * zten _i(ji,jj)1271 zs12(ji,jj) = zfac * z1_ecc2 * pshear_i(ji,jj)1277 zs2(ji,jj) = zfac * z1_ecc2 * ztens(ji,jj) 1278 zs12(ji,jj) = zfac * z1_ecc2 * zshear(ji,jj) * 0.5_wp ! Bug 12 nov 1272 1279 END DO 1273 1280 END DO … … 1347 1354 ! --- Divergence, shear and strength --- ! 1348 1355 IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * zmsk00 ) ! divergence 1349 IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * zmsk00 ) ! shear1356 IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * zmsk00 ) ! maximum shear rate 1350 1357 IF( iom_use('icedlt') ) CALL iom_put( 'icedlt' , pdelta_i * zmsk00 ) ! delta 1351 1358 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength … … 1367 1374 DO ji = 2, jpi - 1 1368 1375 ! Stress invariants 1369 zsig_I(ji,jj) = zs1(ji,jj) * 0.5_wp 1370 zsig_II(ji,jj) = SQRT ( zs2(ji,jj) * zs2(ji,jj) * 0.25_wp + zs12(ji,jj)*zs12(ji,jj) )! 2nd invariant, aka maximum shear stress1376 zsig_I(ji,jj) = zs1(ji,jj) * 0.5_wp ! 1st invariant, aka average normal stress aka negative pressure 1377 zsig_II(ji,jj) = 0.5_wp * SQRT ( zs2(ji,jj) * zs2(ji,jj) + 4. * zs12(ji,jj) * zs12(ji,jj) ) ! 2nd invariant, aka maximum shear stress 1371 1378 END DO 1372 1379 END DO … … 1400 1407 ! following Lemieux & Dupont (2020) 1401 1408 zfac = zp_deltastar_t(ji,jj) 1402 zsig1 = zfac * ( pdivu_i(ji,jj) - zdelta star_t(ji,jj) )1403 zsig2 = zfac * z1_ecc2 * zten _i(ji,jj)1404 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj)1409 zsig1 = zfac * ( pdivu_i(ji,jj) - zdeltat(ji,jj) ) 1410 zsig2 = zfac * z1_ecc2 * ztens(ji,jj) 1411 zsig12 = zfac * z1_ecc2 * zshear(ji,jj) * 0.5_wp ! Bugfix 12 Nov 1405 1412 1406 1413 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 1407 zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st invariant1408 zsig_II(ji,jj) = SQRT ( zsig2 * zsig2 * 0.25_wp +zsig12 * zsig12 ) ! 2nd invariant1414 zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st invariant 1415 zsig_II(ji,jj) = 0.5_wp * SQRT ( zsig2 * zsig2 + 4. *zsig12 * zsig12 ) ! 2nd invariant 1409 1416 1410 1417 ! Normalized principal stresses (used to display the ellipse) … … 1848 1855 ! ! ==================== ! 1849 1856 1850 it_inn_file = ( kt - 1) * kitinntotmax + kitinntot ! time step in the file1851 it_out_file = ( kt - 1 ) * kitoutmax+ kitout1857 it_inn_file = ( kt - nit000 ) * kitinntotmax + kitinntot ! time step in the file 1858 it_out_file = ( kt - nit000 ) * kitoutmax + kitout 1852 1859 1853 1860 ! write variables … … 1887 1894 1888 1895 ! close file 1889 IF( kt == nitend ) istatus = NF90_CLOSE( ncvgid ) 1896 ! IF( kt == nitend ) istatus = NF90_CLOSE( ncvgid ) ! BUG Ed blockley 1897 IF( kt == nitend - nn_fsbc + 1 .AND. kitinntot == kitinntotmax ) istatus = NF90_CLOSE( ncvgid ) 1890 1898 ENDIF 1891 1899
Note: See TracChangeset
for help on using the changeset viewer.
