- Timestamp:
- 2021-06-18T15:21:42+02:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/ticket2680_C1D_PAPA/src/ICE/icedyn_rhg_eap.F90
r14433 r15020 58 58 59 59 REAL(wp), DIMENSION(nx_yield, ny_yield, na_yield) :: s11r, s12r, s22r, s11s, s12s, s22s 60 61 !! * Substitutions 62 # include "do_loop_substitute.h90" 63 # include "domzgr_substitute.h90" 60 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fimask ! mask at F points for the ice 64 61 65 62 !! for convergence tests 66 63 INTEGER :: ncvgid ! netcdf file id 67 64 INTEGER :: nvarid ! netcdf variable id 68 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: aimsk00 69 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: eap_res , aimsk15 65 66 !! * Substitutions 67 # include "do_loop_substitute.h90" 70 68 !!---------------------------------------------------------------------- 71 69 !! NEMO/ICE 4.0 , NEMO Consortium (2018) … … 180 178 REAL(wp), DIMENSION(jpi,jpj) :: ztaux_base, ztauy_base ! ice-bottom stress at U-V points (landfast) 181 179 ! 180 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00, zmsk15 182 181 REAL(wp), DIMENSION(jpi,jpj) :: zmsk01x, zmsk01y ! dummy arrays 183 182 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00x, zmsk00y ! mask for ice presence 184 REAL(wp), DIMENSION(jpi,jpj) :: zfmask ! mask at F points for the ice185 183 186 184 REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter … … 203 201 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_rhg_eap: EAP sea-ice rheology' 204 202 ! 205 IF( kt == nit000 ) THEN 206 ! 207 ! for diagnostics 208 ALLOCATE( aimsk00(jpi,jpj) ) 209 ! for convergence tests 210 IF( nn_rhg_chkcvg > 0 ) ALLOCATE( eap_res(jpi,jpj), aimsk15(jpi,jpj) ) 211 ENDIF 212 ! 203 ! for diagnostics and convergence tests 213 204 DO_2D( 1, 1, 1, 1 ) 214 aimsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice205 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice 215 206 END_2D 216 207 IF( nn_rhg_chkcvg > 0 ) THEN 217 208 DO_2D( 1, 1, 1, 1 ) 218 aimsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less209 zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less 219 210 END_2D 220 211 ENDIF 221 212 ! 222 !!gm for Clem: OPTIMIZATION: I think zfmask can be computed one for all at the initialization....223 213 !------------------------------------------------------------------------------! 224 214 ! 0) mask at F points for the ice 225 215 !------------------------------------------------------------------------------! 226 ! ocean/land mask 227 DO_2D( 1, 0, 1, 0 ) 228 zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) 229 END_2D 230 CALL lbc_lnk( 'icedyn_rhg_eap', zfmask, 'F', 1._wp ) 231 232 ! Lateral boundary conditions on velocity (modify zfmask) 233 DO_2D( 0, 0, 0, 0 ) 234 IF( zfmask(ji,jj) == 0._wp ) THEN 235 zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), & 236 & vmask(ji,jj,1), vmask(ji+1,jj,1) ) ) 216 IF( kt == nit000 ) THEN 217 ! ocean/land mask 218 ALLOCATE( fimask(jpi,jpj) ) 219 IF( rn_ishlat == 0._wp ) THEN 220 DO_2D( 0, 0, 0, 0 ) 221 fimask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) 222 END_2D 223 ELSE 224 DO_2D( 0, 0, 0, 0 ) 225 fimask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) 226 ! Lateral boundary conditions on velocity (modify fimask) 227 IF( fimask(ji,jj) == 0._wp ) THEN 228 fimask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), & 229 & vmask(ji,jj,1), vmask(ji+1,jj,1) ) ) 230 ENDIF 231 END_2D 237 232 ENDIF 238 END_2D 239 DO jj = 2, jpjm1 240 IF( zfmask(1,jj) == 0._wp ) THEN 241 zfmask(1 ,jj) = rn_ishlat * MIN( 1._wp , MAX( vmask(2,jj,1), umask(1,jj+1,1), umask(1,jj,1) ) ) 242 ENDIF 243 IF( zfmask(jpi,jj) == 0._wp ) THEN 244 zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(jpi,jj+1,1), vmask(jpim1,jj,1), umask(jpi,jj-1,1) ) ) 245 ENDIF 246 END DO 247 DO ji = 2, jpim1 248 IF( zfmask(ji,1) == 0._wp ) THEN 249 zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,1,1), umask(ji,2,1), vmask(ji,1,1) ) ) 250 ENDIF 251 IF( zfmask(ji,jpj) == 0._wp ) THEN 252 zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,jpj,1), vmask(ji-1,jpj,1), umask(ji,jpjm1,1) ) ) 253 ENDIF 254 END DO 255 CALL lbc_lnk( 'icedyn_rhg_eap', zfmask, 'F', 1.0_wp ) 233 CALL lbc_lnk( 'icedyn_rhg_eap', fimask, 'F', 1.0_wp ) 234 ENDIF 256 235 257 236 !------------------------------------------------------------------------------! … … 401 380 zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 402 381 & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 403 & ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)382 & ) * r1_e1e2f(ji,jj) * fimask(ji,jj) 404 383 405 384 END_2D … … 760 739 761 740 ! convergence test 762 IF( nn_rhg_chkcvg == 2 ) CALL rhg_cvg_eap( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice )741 IF( nn_rhg_chkcvg == 2 ) CALL rhg_cvg_eap( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice, zmsk15 ) 763 742 ! 764 743 ! ! ==================== ! … … 777 756 zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 778 757 & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 779 & ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)758 & ) * r1_e1e2f(ji,jj) * fimask(ji,jj) 780 759 781 760 END_2D … … 830 809 & ztauy_ai, 'V', -1.0_wp, ztaux_bi, 'U', -1.0_wp, ztauy_bi, 'V', -1.0_wp ) 831 810 ! 832 CALL iom_put( 'utau_oi' , ztaux_oi * aimsk00 )833 CALL iom_put( 'vtau_oi' , ztauy_oi * aimsk00 )834 CALL iom_put( 'utau_ai' , ztaux_ai * aimsk00 )835 CALL iom_put( 'vtau_ai' , ztauy_ai * aimsk00 )836 CALL iom_put( 'utau_bi' , ztaux_bi * aimsk00 )837 CALL iom_put( 'vtau_bi' , ztauy_bi * aimsk00 )811 CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 ) 812 CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 ) 813 CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 ) 814 CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 ) 815 CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 ) 816 CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 ) 838 817 ENDIF 839 818 840 819 ! --- divergence, shear and strength --- ! 841 IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * aimsk00 ) ! divergence842 IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * aimsk00 ) ! shear843 IF( iom_use('icedlt') ) CALL iom_put( 'icedlt' , pdelta_i * aimsk00 ) ! delta844 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * aimsk00 ) ! strength820 IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * zmsk00 ) ! divergence 821 IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * zmsk00 ) ! shear 822 IF( iom_use('icedlt') ) CALL iom_put( 'icedlt' , pdelta_i * zmsk00 ) ! delta 823 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength 845 824 846 825 ! --- Stress tensor invariants (SIMIP diags) --- ! … … 867 846 ! 868 847 ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008) 869 IF( iom_use('normstr') ) CALL iom_put( 'normstr', zsig_I (:,:) * aimsk00(:,:) ) ! Normal stress870 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr', zsig_II(:,:) * aimsk00(:,:) ) ! Maximum shear stress848 IF( iom_use('normstr') ) CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress 849 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress 871 850 872 851 DEALLOCATE ( zsig_I, zsig_II ) … … 914 893 CALL lbc_lnk( 'icedyn_rhg_eap', zyield11, 'T', 1.0_wp, zyield22, 'T', 1.0_wp, zyield12, 'T', 1.0_wp ) 915 894 916 CALL iom_put( 'yield11', zyield11 * aimsk00 )917 CALL iom_put( 'yield22', zyield22 * aimsk00 )918 CALL iom_put( 'yield12', zyield12 * aimsk00 )895 CALL iom_put( 'yield11', zyield11 * zmsk00 ) 896 CALL iom_put( 'yield22', zyield22 * zmsk00 ) 897 CALL iom_put( 'yield12', zyield12 * zmsk00 ) 919 898 ENDIF 920 899 … … 922 901 IF( iom_use('aniso') ) THEN 923 902 CALL lbc_lnk( 'icedyn_rhg_eap', paniso_11, 'T', 1.0_wp ) 924 CALL iom_put( 'aniso' , paniso_11 * aimsk00 )903 CALL iom_put( 'aniso' , paniso_11 * zmsk00 ) 925 904 ENDIF 926 905 … … 933 912 & zfU, 'U', -1.0_wp, zfV, 'V', -1.0_wp ) 934 913 935 CALL iom_put( 'dssh_dx' , zspgU * aimsk00 ) ! Sea-surface tilt term in force balance (x)936 CALL iom_put( 'dssh_dy' , zspgV * aimsk00 ) ! Sea-surface tilt term in force balance (y)937 CALL iom_put( 'corstrx' , zCorU * aimsk00 ) ! Coriolis force term in force balance (x)938 CALL iom_put( 'corstry' , zCorV * aimsk00 ) ! Coriolis force term in force balance (y)939 CALL iom_put( 'intstrx' , zfU * aimsk00 ) ! Internal force term in force balance (x)940 CALL iom_put( 'intstry' , zfV * aimsk00 ) ! Internal force term in force balance (y)914 CALL iom_put( 'dssh_dx' , zspgU * zmsk00 ) ! Sea-surface tilt term in force balance (x) 915 CALL iom_put( 'dssh_dy' , zspgV * zmsk00 ) ! Sea-surface tilt term in force balance (y) 916 CALL iom_put( 'corstrx' , zCorU * zmsk00 ) ! Coriolis force term in force balance (x) 917 CALL iom_put( 'corstry' , zCorV * zmsk00 ) ! Coriolis force term in force balance (y) 918 CALL iom_put( 'intstrx' , zfU * zmsk00 ) ! Internal force term in force balance (x) 919 CALL iom_put( 'intstry' , zfV * zmsk00 ) ! Internal force term in force balance (y) 941 920 ENDIF 942 921 … … 949 928 DO_2D( 0, 0, 0, 0 ) 950 929 ! 2D ice mass, snow mass, area transport arrays (X, Y) 951 zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * aimsk00(ji,jj)952 zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * aimsk00(ji,jj)930 zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj) 931 zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj) 953 932 954 933 zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component … … 984 963 IF( ln_aEVP ) THEN ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 985 964 CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , & 986 & ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * aimsk15(:,:) )965 & ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) ) 987 966 ELSE ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 988 967 CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , & 989 & ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * aimsk15(:,:) )968 & ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) ) 990 969 ENDIF 991 970 ENDIF … … 995 974 996 975 997 SUBROUTINE rhg_cvg_eap( kt, kiter, kitermax, pu, pv, pub, pvb )976 SUBROUTINE rhg_cvg_eap( kt, kiter, kitermax, pu, pv, pub, pvb, pmsk15 ) 998 977 !!---------------------------------------------------------------------- 999 978 !! *** ROUTINE rhg_cvg_eap *** … … 1010 989 INTEGER , INTENT(in) :: kt, kiter, kitermax ! ocean time-step index 1011 990 REAL(wp), DIMENSION(:,:), INTENT(in) :: pu, pv, pub, pvb ! now and before velocities 991 REAL(wp), DIMENSION(:,:), INTENT(in) :: pmsk15 1012 992 !! 1013 993 INTEGER :: it, idtime, istatus … … 1044 1024 zresm = 0._wp 1045 1025 ELSE 1046 DO_2D( 1, 1, 1, 1 ) 1047 eap_res(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 1048 & ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * aimsk15(ji,jj) 1026 zresm = 0._wp 1027 DO_2D( 0, 0, 0, 0 ) 1028 zresm = MAX( zresm, MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 1029 & ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * pmsk15(ji,jj) ) 1049 1030 END_2D 1050 1051 zresm = MAXVAL( eap_res )1052 1031 CALL mpp_max( 'icedyn_rhg_evp', zresm ) ! max over the global domain 1053 1032 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
