Changeset 13799
- Timestamp:
- 2020-11-16T12:04:05+01:00 (4 years ago)
- Location:
- NEMO/branches/2020/dev_r13787_SI3-10_EAP
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13787_SI3-10_EAP/src/ICE/icedyn_rhg_eap.F90
r13797 r13799 430 430 & + zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) & 431 431 & ) * 0.25_wp * r1_e1e2t(ji,jj) 432 !WRITE(numout,*) 'shear output', ji, jj, zdsT433 432 434 433 ! divergence at T points (duplication to avoid communications) … … 447 446 448 447 ! structure tensor update 449 !!$ IF (mod(jter,10) == 0) THEN450 448 CALL calc_ffrac(zstressptmp, zstressmtmp, zstress12tmp(ji,jj), paniso_11(ji,jj), paniso_12(ji,jj), zmresult11, zmresult12) 451 449 452 !!$ paniso_11(ji,jj) = (paniso_11(ji,jj)*z1_dtevp_A + zp5kth + zmresult11) * z1dtevpkth ! implicit453 !!$ paniso_12(ji,jj) = (paniso_12(ji,jj)*z1_dtevp_A + zmresult12) * z1dtevpkth ! implicit454 450 paniso_11(ji,jj) = (paniso_11(ji,jj) + 0.5*2.e-5*zdtevp + zmresult11*zdtevp) / (1. + 2.e-5*zdtevp) ! implicit 455 451 paniso_12(ji,jj) = (paniso_12(ji,jj) + zmresult12*zdtevp) / (1. + 2.e-5*zdtevp) ! implicit 456 !!$ ENDIF457 452 458 453 IF (jter == nn_nevp) THEN … … 479 474 480 475 ! stress at T points (zkt/=0 if landfast) 481 !zs1(ji,jj) = ( zs1(ji,jj)*zalph1 + zp_delt(ji,jj) * ( zdiv*(1._wp + zkt) - zdelta(ji,jj)*(1._wp - zkt) ) ) * z1_alph1482 !zs2(ji,jj) = ( zs2(ji,jj)*zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2483 !!$ zs1(ji,jj) = ( zs1(ji,jj) + zstressptmp * zalph1 ) * z1_alph1484 !!$ zs2(ji,jj) = ( zs2(ji,jj) + zstressmtmp * zalph1 ) * z1_alph1485 476 zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zstressptmp ) * z1_alph1 486 477 zs2(ji,jj) = ( zs2(ji,jj) * zalph1 + zstressmtmp ) * z1_alph1 487 !zs1(ji,jj) = ( stressptmp * zs1(ji,jj) + zalph1 ) * z1_alph1488 !zs2(ji,jj) = ( stressmtmp * zs2(ji,jj) + zalph1 ) * z1_alph1489 !WRITE(numout,*) 'stress output', ji, jj, zs1(ji,jj)490 478 END_2D 491 479 CALL lbc_lnk_multi( 'icedyn_rhg_eap', zstress12tmp, 'T', 1.0_wp , paniso_11, 'T', 1.0_wp , paniso_12, 'T', 1.0_wp) … … 514 502 515 503 ! stress at F points (zkt/=0 if landfast) 516 !zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2517 !!$ zs12(ji,jj) = ( zs12(ji,jj) + zstress12tmpF * zalph1 ) * z1_alph1518 504 zs12(ji,jj) = ( zs12(ji,jj) * zalph1 + zstress12tmpF ) * z1_alph1 519 !zs12(ji,jj) = ( stress12tmpF * zs12(ji,jj) + zalph1 ) * z1_alph1520 505 521 506 END_2D … … 598 583 ! 599 584 #if defined key_agrif 600 !! 585 !! CALL agrif_interp_ice( 'V', jter, nn_nevp ) 601 586 CALL agrif_interp_ice( 'V' ) 602 587 #endif … … 650 635 ! 651 636 #if defined key_agrif 652 !! 637 !! CALL agrif_interp_ice( 'U', jter, nn_nevp ) 653 638 CALL agrif_interp_ice( 'U' ) 654 639 #endif … … 704 689 ! 705 690 #if defined key_agrif 706 !! 691 !! CALL agrif_interp_ice( 'U', jter, nn_nevp ) 707 692 CALL agrif_interp_ice( 'U' ) 708 693 #endif … … 756 741 ! 757 742 #if defined key_agrif 758 !! 743 !! CALL agrif_interp_ice( 'V', jter, nn_nevp ) 759 744 CALL agrif_interp_ice( 'V' ) 760 745 #endif … … 1277 1262 zsig22 = pstrength*(zstemp22r + kfriction*zstemp22s) * zinvsin 1278 1263 1279 !WRITE(numout,*) 'principal axis stress inside loop', sig11, sig12, sig221280 1281 1264 ! Back - rotation of the stress from principal axes into general coordinates 1282 1265 … … 1438 1421 ! 1439 1422 IF( MIN( id1, id2, id3, id4, id5 ) > 0 ) THEN ! fields exist 1440 CALL iom_get( numrir, jpdom_auto glo, 'stress1_i' , stress1_i)1441 CALL iom_get( numrir, jpdom_auto glo, 'stress2_i' , stress2_i)1442 CALL iom_get( numrir, jpdom_auto glo, 'stress12_i', stress12_i)1443 CALL iom_get( numrir, jpdom_auto glo, 'aniso_11' , aniso_11)1444 CALL iom_get( numrir, jpdom_auto glo, 'aniso_12' , aniso_12)1423 CALL iom_get( numrir, jpdom_auto, 'stress1_i' , stress1_i , cd_type = 'T' ) 1424 CALL iom_get( numrir, jpdom_auto, 'stress2_i' , stress2_i , cd_type = 'T' ) 1425 CALL iom_get( numrir, jpdom_auto, 'stress12_i', stress12_i, cd_type = 'F' ) 1426 CALL iom_get( numrir, jpdom_auto, 'aniso_11' , aniso_11 , cd_type = 'T' ) 1427 CALL iom_get( numrir, jpdom_auto, 'aniso_12' , aniso_12 , cd_type = 'T' ) 1445 1428 ELSE ! start rheology from rest 1446 1429 IF(lwp) WRITE(numout,*) -
NEMO/branches/2020/dev_r13787_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/icedyn_rhg_eap.F90
r13797 r13799 13 13 !! 3.7 ! 2017 (C. Rousset) add aEVP (Kimmritz 2016-2017) 14 14 !! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube] 15 !! ! 2019 (S. Rynders ) change into eap rheology from15 !! ! 2019 (S. Rynders, Y. Aksenov, C. Rousset) change into eap rheology from 16 16 !! CICE code (Tsamados, Heorton) 17 17 !!---------------------------------------------------------------------- … … 450 450 451 451 ! structure tensor update 452 !!$ IF (mod(jter,10) == 0) THEN453 452 CALL calc_ffrac(zstressptmp, zstressmtmp, zstress12tmp(ji,jj), paniso_11(ji,jj), paniso_12(ji,jj), zmresult11, zmresult12) 454 453 455 !!$ paniso_11(ji,jj) = (paniso_11(ji,jj)*z1_dtevp_A + zp5kth + zmresult11) * z1dtevpkth ! implicit456 !!$ paniso_12(ji,jj) = (paniso_12(ji,jj)*z1_dtevp_A + zmresult12) * z1dtevpkth ! implicit457 454 paniso_11(ji,jj) = (paniso_11(ji,jj) + 0.5*2.e-5*zdtevp + zmresult11*zdtevp) / (1. + 2.e-5*zdtevp) ! implicit 458 455 paniso_12(ji,jj) = (paniso_12(ji,jj) + zmresult12*zdtevp) / (1. + 2.e-5*zdtevp) ! implicit 459 !!$ ENDIF460 456 461 457 IF (jter == nn_nevp) THEN … … 510 506 511 507 ! stress at F points (zkt/=0 if landfast) 512 !zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2513 !!$ zs12(ji,jj) = ( zs12(ji,jj) + zstress12tmpF * zalph1 ) * z1_alph1514 508 zs12(ji,jj) = ( zs12(ji,jj) * zalph1 + zstress12tmpF ) * z1_alph1 515 !zs12(ji,jj) = ( stress12tmpF * zs12(ji,jj) + zalph1 ) * z1_alph1516 509 517 510 END_2D … … 599 592 ! 600 593 #if defined key_agrif 601 !! 594 !! CALL agrif_interp_ice( 'V', jter, nn_nevp ) 602 595 CALL agrif_interp_ice( 'V' ) 603 596 #endif … … 656 649 ! 657 650 #if defined key_agrif 658 !! 651 !! CALL agrif_interp_ice( 'U', jter, nn_nevp ) 659 652 CALL agrif_interp_ice( 'U' ) 660 653 #endif … … 714 707 ! 715 708 #if defined key_agrif 716 !! 709 !! CALL agrif_interp_ice( 'U', jter, nn_nevp ) 717 710 CALL agrif_interp_ice( 'U' ) 718 711 #endif … … 770 763 ! 771 764 #if defined key_agrif 772 !! 765 !! CALL agrif_interp_ice( 'V', jter, nn_nevp ) 773 766 CALL agrif_interp_ice( 'V' ) 774 767 #endif … … 1458 1451 CALL iom_get( numrir, jpdom_auto, 'stress2_i' , stress2_i , cd_type = 'T' ) 1459 1452 CALL iom_get( numrir, jpdom_auto, 'stress12_i', stress12_i, cd_type = 'F' ) 1460 CALL iom_get( numrir, jpdom_auto, 'aniso_11' , aniso_11 , cd_type = 'T' 1461 CALL iom_get( numrir, jpdom_auto, 'aniso_12' , aniso_12 , cd_type = 'T' 1453 CALL iom_get( numrir, jpdom_auto, 'aniso_11' , aniso_11 , cd_type = 'T' ) 1454 CALL iom_get( numrir, jpdom_auto, 'aniso_12' , aniso_12 , cd_type = 'T' ) 1462 1455 ELSE ! start rheology from rest 1463 1456 IF(lwp) WRITE(numout,*) -
NEMO/branches/2020/dev_r13787_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/usrdef_hgr.F90
r13797 r13799 28 28 !! * Substitutions 29 29 # include "do_loop_substitute.h90" 30 # include "domzgr_substitute.h90"31 30 !!---------------------------------------------------------------------- 32 31 !! NEMO/OCE 4.0 , NEMO Consortium (2018) -
NEMO/branches/2020/dev_r13787_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/usrdef_sbc.F90
r13797 r13799 34 34 35 35 !! * Substitutions 36 # include "do_loop_substitute.h90" 36 37 !!---------------------------------------------------------------------- 37 38 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 119 120 IF( kt==nit000 .AND. lwp) WRITE(numout,*)' usrdef_sbc_ice : ICE_RHEO case: analytical stress forcing' 120 121 121 DO jj = 2, jpjm1 122 DO ji = 2, jpim1 123 ! wind spins up over 6 hours, factor 1000 to balance the units 124 windu(ji,jj) = Umax/sqrt(d*1000)*(d-2*mig(ji)*res)/((d-2*mig(ji)*res)**2+(d-2*mjg(jj)*res)**2*Rwind**2)**(1/4)*min(kt*30./21600,1.) 125 windv(ji,jj) = Umax/sqrt(d*1000)*(d-2*mjg(jj)*res)/((d-2*mig(ji)*res)**2+(d-2*mjg(jj)*res)**2*Rwind**2)**(1/4)*Rwind*min(kt*30./21600,1.) 126 END DO 127 END DO 122 DO_2D( 0, 0, 0, 0 ) 123 ! wind spins up over 6 hours, factor 1000 to balance the units 124 windu(ji,jj) = Umax/sqrt(d*1000)*(d-2*mig(ji)*res)/((d-2*mig(ji)*res)**2+(d-2*mjg(jj)*res)**2*Rwind**2)**(1/4)*min(kt*30./21600,1.) 125 windv(ji,jj) = Umax/sqrt(d*1000)*(d-2*mjg(jj)*res)/((d-2*mig(ji)*res)**2+(d-2*mjg(jj)*res)**2*Rwind**2)**(1/4)*Rwind*min(kt*30./21600,1.) 126 END_2D 128 127 CALL lbc_lnk_multi( 'usrdef_sbc', windu, 'U', -1., windv, 'V', -1. ) 129 128 … … 134 133 ! ------------------------------------------------------------ ! 135 134 ! C-grid ice dynamics : U & V-points (same as ocean) 136 DO jj = 2, jpjm1 137 DO ji = 2, jpim1 138 zwndi_t = ( windu(ji,jj) - r_vfac * 0.5 * ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) ) 139 zwndj_t = ( windv(ji,jj) - r_vfac * 0.5 * ( v_ice(ji,jj-1) + v_ice(ji,jj) ) ) 140 wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 141 END DO 142 END DO 135 DO_2D( 0, 0, 0, 0 ) 136 zwndi_t = ( windu(ji,jj) - r_vfac * 0.5 * ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) ) 137 zwndj_t = ( windv(ji,jj) - r_vfac * 0.5 * ( v_ice(ji,jj-1) + v_ice(ji,jj) ) ) 138 wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 139 END_2D 143 140 CALL lbc_lnk( 'usrdef_sbc', wndm_ice, 'T', 1. ) 144 141 … … 152 149 ! ------------------------------------------------------------ ! 153 150 ! C-grid ice dynamics : U & V-points (same as ocean) 154 DO jj = 2, jpjm1 155 DO ji = 2, jpim1 156 utau_ice(ji,jj) = 0.5 * zrhoa * Cd_atm * ( wndm_ice(ji+1,jj ) + wndm_ice(ji,jj) ) & 157 & * ( 0.5 * (windu(ji+1,jj) + windu(ji,jj) ) - r_vfac * u_ice(ji,jj) ) 158 vtau_ice(ji,jj) = 0.5 * zrhoa * Cd_atm * ( wndm_ice(ji,jj+1 ) + wndm_ice(ji,jj) ) & 159 & * ( 0.5 * (windv(ji,jj+1) + windv(ji,jj) ) - r_vfac * v_ice(ji,jj) ) 160 END DO 161 END DO 151 DO_2D( 0, 0, 0, 0 ) 152 utau_ice(ji,jj) = 0.5 * zrhoa * Cd_atm * ( wndm_ice(ji+1,jj ) + wndm_ice(ji,jj) ) & 153 & * ( 0.5 * (windu(ji+1,jj) + windu(ji,jj) ) - r_vfac * u_ice(ji,jj) ) 154 vtau_ice(ji,jj) = 0.5 * zrhoa * Cd_atm * ( wndm_ice(ji,jj+1 ) + wndm_ice(ji,jj) ) & 155 & * ( 0.5 * (windv(ji,jj+1) + windv(ji,jj) ) - r_vfac * v_ice(ji,jj) ) 156 END_2D 162 157 CALL lbc_lnk_multi( 'usrdef_sbc', utau_ice, 'U', -1., vtau_ice, 'V', -1. ) 163 158 !
Note: See TracChangeset
for help on using the changeset viewer.