- Timestamp:
- 2017-09-05T19:53:41+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90
r8486 r8498 31 31 USE prtctl ! Print control 32 32 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 33 USE iom ! 33 34 #if defined key_agrif 34 35 USE agrif_lim3_interp … … 111 112 INTEGER :: ji, jj ! dummy loop indices 112 113 INTEGER :: jter ! local integers 113 CHARACTER (len=50) :: charout114 114 115 115 REAL(wp) :: zrhoco ! rau0 * rn_cio … … 121 121 REAL(wp) :: zTauO, zTauB, zTauE, zvel ! temporary scalars 122 122 123 REAL(wp) :: zsig1, zsig2 ! internal ice stress124 123 REAL(wp) :: zresm ! Maximal error on ice velocity 125 124 REAL(wp) :: zintb, zintn ! dummy argument 126 125 REAL(wp) :: zfac_x, zfac_y 126 REAL(wp) :: zshear, zdum1, zdum2 127 127 128 128 REAL(wp), DIMENSION(jpi,jpj) :: z1_e1t0, z1_e2t0 ! scale factors … … 152 152 REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter 153 153 REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity equals ocean velocity 154 !! --- diags 155 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zswi, zsig1, zsig2, zsig3 156 !! --- SIMIP diags 157 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_sig1 ! Average normal stress in sea ice 158 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_sig2 ! Maximum shear stress in sea ice 159 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_dssh_dx ! X-direction sea-surface tilt term (N/m2) 160 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_dssh_dy ! X-direction sea-surface tilt term (N/m2) 161 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_corstrx ! X-direction coriolis stress (N/m2) 162 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_corstry ! Y-direction coriolis stress (N/m2) 163 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_intstrx ! X-direction internal stress (N/m2) 164 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_intstry ! Y-direction internal stress (N/m2) 165 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_utau_oi ! X-direction ocean-ice stress 166 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_vtau_oi ! Y-direction ocean-ice stress 167 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) 168 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s) 169 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s) 170 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s) 171 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xatrp ! X-component of area transport (m2/s) 172 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_yatrp ! Y-component of area transport (m2/s) 154 173 !!------------------------------------------------------------------- 155 174 … … 671 690 672 691 !------------------------------------------------------------------------------! 673 ! 5) SIMIP diagnostics 674 !------------------------------------------------------------------------------! 675 676 !!gm encapsulate with a flag (iom_use of the variable or better a flag defined one for all testing if one of the 677 !! diag is output. NB the diag_... are should only be ALLOCATED if the flag is true ! 678 679 DO jj = 2, jpjm1 680 DO ji = 2, jpim1 681 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 682 683 ! Stress tensor invariants (normal and shear stress N/m) 684 diag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch ! normal stress 685 diag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch ! shear stress 686 687 ! Stress terms of the momentum equation (N/m2) 688 diag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch ! sea surface slope stress term 689 diag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch 690 691 diag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch ! Coriolis stress term 692 diag_corstry(ji,jj) = zCory(ji,jj) * rswitch 693 694 diag_intstrx(ji,jj) = zfU(ji,jj) * rswitch ! internal stress term 695 diag_intstry(ji,jj) = zfV(ji,jj) * rswitch 696 697 diag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch ! oceanic stress 698 diag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch 699 700 ! 2D ice mass, snow mass, area transport arrays (X, Y) 701 zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch 702 zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch 703 704 diag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component 705 diag_ymtrp_ice(ji,jj) = rhoic * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) ! '' Y- '' 706 707 diag_xmtrp_snw(ji,jj) = rhosn * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component 708 diag_ymtrp_snw(ji,jj) = rhosn * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) ! '' Y- '' 709 710 diag_xatrp(ji,jj) = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) ) ! area transport, X-component 711 diag_yatrp(ji,jj) = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) ) ! '' Y- '' 712 713 END DO 714 END DO 715 716 CALL lbc_lnk_multi( diag_sig1 , 'T', 1., diag_sig2 , 'T', 1., & 717 & diag_dssh_dx, 'U', -1., diag_dssh_dy, 'V', -1., & 718 & diag_corstrx, 'U', -1., diag_corstry, 'V', -1., & 719 & diag_intstrx, 'U', -1., diag_intstry, 'V', -1. ) 720 721 CALL lbc_lnk_multi( diag_utau_oi, 'U', -1., diag_vtau_oi, 'V', -1. ) 722 723 CALL lbc_lnk_multi( diag_xmtrp_ice, 'U', -1., diag_xmtrp_snw, 'U', -1., & 724 & diag_xatrp , 'U', -1., diag_ymtrp_ice, 'V', -1., & 725 & diag_ymtrp_snw, 'V', -1., diag_yatrp , 'V', -1. ) 726 727 ! 728 !------------------------------------------------------------------------------! 729 ! 6) Control prints of residual and charge ellipse 730 !------------------------------------------------------------------------------! 731 ! 732 ! print the residual for convergence 733 IF(ln_ctl) THEN 734 WRITE(charout,FMT="('ice_rhg_evp : res =',D23.16, ' iter =',I4)") zresm, jter 735 CALL prt_ctl_info(charout) 736 CALL prt_ctl(tab2d_1=u_ice, clinfo1=' ice_rhg_evp : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :') 692 ! 5) diagnostics 693 !------------------------------------------------------------------------------! 694 ! --- ellipse --- ! 695 IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') ) THEN 737 696 ! 738 ! print charge ellipse 739 ! This can be desactivated once the user is sure that the stress state 740 ! lie on the charge ellipse. See Bouillon et al. (2008) for more details 741 IF( MOD(kt_ice+nn_fsbc-1,nwrite) == 0 ) THEN 742 WRITE(charout,FMT="('ice_rhg_evp :', I4, I6, I1, I1, A10)") 1000, kt_ice, 0, 0, ' ch. ell. ' 743 CALL prt_ctl_info(charout) 744 DO jj = 2, jpjm1 745 DO ji = 2, jpim1 746 IF (strength(ji,jj) > 1.0) THEN 747 zsig1 = ( zs1(ji,jj) + SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*strength(ji,jj) ) 748 zsig2 = ( zs1(ji,jj) - SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*strength(ji,jj) ) 749 WRITE(charout,FMT="('ice_rhg_evp :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)") 750 CALL prt_ctl_info(charout) 751 ENDIF 752 END DO 753 END DO 754 WRITE(charout,FMT="('ice_rhg_evp :', I4, I6, I1, I1, A10)") 2000, kt_ice, 0, 0, ' ch. ell. ' 755 CALL prt_ctl_info(charout) 756 ENDIF 697 ALLOCATE( zswi(jpi,jpj) , zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) ) 698 ! 699 DO jj = 1, jpj 700 DO ji = 1, jpi 701 zswi(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 702 END DO 703 END DO 704 705 DO jj = 2, jpjm1 706 DO ji = 2, jpim1 707 zdum1 = ( zswi(ji-1,jj) * stress12_i(ji-1,jj) + zswi(ji ,jj-1) * stress12_i(ji ,jj-1) + & ! stress12_i at T-point 708 & zswi(ji ,jj) * stress12_i(ji ,jj) + zswi(ji-1,jj-1) * stress12_i(ji-1,jj-1) ) & 709 & / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) ) 710 711 zshear = SQRT( stress2_i(ji,jj) * stress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 712 713 zdum2 = zswi(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 714 715 !! zsig1(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 716 !! zsig2(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002) 717 !! zsig3(ji,jj) = zdum2**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress 718 !! ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 719 zsig1(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) ) ! compressive stress, see Bouillon et al. 2015 720 zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear ) ! shear stress 721 zsig3(ji,jj) = zdum2**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 722 END DO 723 END DO 724 CALL lbc_lnk_multi( zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. ) 725 ! 726 IF( iom_use('isig1') ) CALL iom_put( "isig1" , zsig1 ) 727 IF( iom_use('isig2') ) CALL iom_put( "isig2" , zsig2 ) 728 IF( iom_use('isig3') ) CALL iom_put( "isig3" , zsig3 ) 729 ! 730 DEALLOCATE( zswi , zsig1 , zsig2 , zsig3 ) 731 ENDIF 732 733 ! --- SIMIP --- ! 734 IF ( iom_use( 'normstr' ) .OR. iom_use( 'sheastr' ) .OR. iom_use( 'dssh_dx' ) .OR. iom_use( 'dssh_dy' ) .OR. & 735 & iom_use( 'corstrx' ) .OR. iom_use( 'corstry' ) .OR. iom_use( 'intstrx' ) .OR. iom_use( 'intstry' ) .OR. & 736 & iom_use( 'utau_oi' ) .OR. iom_use( 'vtau_oi' ) .OR. iom_use( 'xmtrpice' ) .OR. iom_use( 'ymtrpice' ) .OR. & 737 & iom_use( 'xmtrpsnw' ) .OR. iom_use( 'ymtrpsnw' ) .OR. iom_use( 'xatrp' ) .OR. iom_use( 'yatrp' ) ) THEN 738 739 ALLOCATE( zdiag_sig1 (jpi,jpj) , zdiag_sig2 (jpi,jpj) , zdiag_dssh_dx (jpi,jpj) , zdiag_dssh_dy (jpi,jpj) , & 740 & zdiag_corstrx (jpi,jpj) , zdiag_corstry (jpi,jpj) , zdiag_intstrx (jpi,jpj) , zdiag_intstry (jpi,jpj) , & 741 & zdiag_utau_oi (jpi,jpj) , zdiag_vtau_oi (jpi,jpj) , zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , & 742 & zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp (jpi,jpj) , zdiag_yatrp (jpi,jpj) ) 743 744 DO jj = 2, jpjm1 745 DO ji = 2, jpim1 746 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 747 748 ! Stress tensor invariants (normal and shear stress N/m) 749 zdiag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch ! normal stress 750 zdiag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch ! shear stress 751 752 ! Stress terms of the momentum equation (N/m2) 753 zdiag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch ! sea surface slope stress term 754 zdiag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch 755 756 zdiag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch ! Coriolis stress term 757 zdiag_corstry(ji,jj) = zCory(ji,jj) * rswitch 758 759 zdiag_intstrx(ji,jj) = zfU(ji,jj) * rswitch ! internal stress term 760 zdiag_intstry(ji,jj) = zfV(ji,jj) * rswitch 761 762 zdiag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch ! oceanic stress 763 zdiag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch 764 765 ! 2D ice mass, snow mass, area transport arrays (X, Y) 766 zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch 767 zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch 768 769 zdiag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component 770 zdiag_ymtrp_ice(ji,jj) = rhoic * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) ! '' Y- '' 771 772 zdiag_xmtrp_snw(ji,jj) = rhosn * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component 773 zdiag_ymtrp_snw(ji,jj) = rhosn * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) ! '' Y- '' 774 775 zdiag_xatrp(ji,jj) = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) ) ! area transport, X-component 776 zdiag_yatrp(ji,jj) = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) ) ! '' Y- '' 777 778 END DO 779 END DO 780 781 CALL lbc_lnk_multi( zdiag_sig1 , 'T', 1., zdiag_sig2 , 'T', 1., & 782 & zdiag_dssh_dx, 'U', -1., zdiag_dssh_dy, 'V', -1., & 783 & zdiag_corstrx, 'U', -1., zdiag_corstry, 'V', -1., & 784 & zdiag_intstrx, 'U', -1., zdiag_intstry, 'V', -1. ) 785 786 CALL lbc_lnk_multi( zdiag_utau_oi, 'U', -1., zdiag_vtau_oi, 'V', -1. ) 787 788 CALL lbc_lnk_multi( zdiag_xmtrp_ice, 'U', -1., zdiag_xmtrp_snw, 'U', -1., & 789 & zdiag_xatrp , 'U', -1., zdiag_ymtrp_ice, 'V', -1., & 790 & zdiag_ymtrp_snw, 'V', -1., zdiag_yatrp , 'V', -1. ) 791 792 IF ( iom_use( 'normstr' ) ) CALL iom_put( 'normstr' , zdiag_sig1(:,:) ) ! Normal stress 793 IF ( iom_use( 'sheastr' ) ) CALL iom_put( 'sheastr' , zdiag_sig2(:,:) ) ! Shear stress 794 IF ( iom_use( 'dssh_dx' ) ) CALL iom_put( 'dssh_dx' , zdiag_dssh_dx(:,:) ) ! Sea-surface tilt term in force balance (x) 795 IF ( iom_use( 'dssh_dy' ) ) CALL iom_put( 'dssh_dy' , zdiag_dssh_dy(:,:) ) ! Sea-surface tilt term in force balance (y) 796 IF ( iom_use( 'corstrx' ) ) CALL iom_put( 'corstrx' , zdiag_corstrx(:,:) ) ! Coriolis force term in force balance (x) 797 IF ( iom_use( 'corstry' ) ) CALL iom_put( 'corstry' , zdiag_corstry(:,:) ) ! Coriolis force term in force balance (y) 798 IF ( iom_use( 'intstrx' ) ) CALL iom_put( 'intstrx' , zdiag_intstrx(:,:) ) ! Internal force term in force balance (x) 799 IF ( iom_use( 'intstry' ) ) CALL iom_put( 'intstry' , zdiag_intstry(:,:) ) ! Internal force term in force balance (y) 800 IF ( iom_use( 'utau_oi' ) ) CALL iom_put( 'utau_oi' , zdiag_utau_oi(:,:) ) ! Ocean stress term in force balance (x) 801 IF ( iom_use( 'vtau_oi' ) ) CALL iom_put( 'vtau_oi' , zdiag_vtau_oi(:,:) ) ! Ocean stress term in force balance (y) 802 IF ( iom_use( 'xmtrpice' ) ) CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice(:,:) ) ! X-component of sea-ice mass transport (kg/s) 803 IF ( iom_use( 'ymtrpice' ) ) CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice(:,:) ) ! Y-component of sea-ice mass transport 804 IF ( iom_use( 'xmtrpsnw' ) ) CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw(:,:) ) ! X-component of snow mass transport (kg/s) 805 IF ( iom_use( 'ymtrpsnw' ) ) CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw(:,:) ) ! Y-component of snow mass transport 806 IF ( iom_use( 'xatrp' ) ) CALL iom_put( 'xatrp' , zdiag_xatrp(:,:) ) ! X-component of ice area transport 807 IF ( iom_use( 'yatrp' ) ) CALL iom_put( 'yatrp' , zdiag_yatrp(:,:) ) ! Y-component of ice area transport 808 809 DEALLOCATE( zdiag_sig1 , zdiag_sig2 , zdiag_dssh_dx , zdiag_dssh_dy , & 810 & zdiag_corstrx , zdiag_corstry , zdiag_intstrx , zdiag_intstry , & 811 & zdiag_utau_oi , zdiag_vtau_oi , zdiag_xmtrp_ice , zdiag_ymtrp_ice , & 812 & zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp , zdiag_yatrp ) 813 757 814 ENDIF 758 815 !
Note: See TracChangeset
for help on using the changeset viewer.
