- Timestamp:
- 2021-07-16T20:00:12+02:00 (3 years ago)
- Location:
- NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/ICE/icedyn_rhg_eap.F90 (modified) (16 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO
- Property svn:externals
-
old new 9 9 10 10 # SETTE 11 ^/utils/CI/sette@14244 sette 11 ^/utils/CI/sette@HEAD sette 12
-
- Property svn:externals
-
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/ICE/icedyn_rhg_eap.F90
r14448 r15127 58 58 59 59 REAL(wp), DIMENSION(nx_yield, ny_yield, na_yield) :: s11r, s12r, s22r, s11s, s12s, s22s 60 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fimask ! mask at F points for the ice 61 62 !! for convergence tests 63 INTEGER :: ncvgid ! netcdf file id 64 INTEGER :: nvarid ! netcdf variable id 60 65 61 66 !! * Substitutions 62 67 # include "do_loop_substitute.h90" 63 68 # include "domzgr_substitute.h90" 64 65 !! for convergence tests66 INTEGER :: ncvgid ! netcdf file id67 INTEGER :: nvarid ! netcdf variable id68 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: aimsk0069 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: eap_res , aimsk1570 69 !!---------------------------------------------------------------------- 71 70 !! NEMO/ICE 4.0 , NEMO Consortium (2018) … … 180 179 REAL(wp), DIMENSION(jpi,jpj) :: ztaux_base, ztauy_base ! ice-bottom stress at U-V points (landfast) 181 180 ! 181 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00, zmsk15 182 182 REAL(wp), DIMENSION(jpi,jpj) :: zmsk01x, zmsk01y ! dummy arrays 183 183 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 184 186 185 REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter … … 203 202 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_rhg_eap: EAP sea-ice rheology' 204 203 ! 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 ! 204 ! for diagnostics and convergence tests 213 205 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 ice206 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice 215 207 END_2D 216 208 IF( nn_rhg_chkcvg > 0 ) THEN 217 209 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 less210 zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less 219 211 END_2D 220 212 ENDIF 221 213 ! 222 !!gm for Clem: OPTIMIZATION: I think zfmask can be computed one for all at the initialization....223 214 !------------------------------------------------------------------------------! 224 215 ! 0) mask at F points for the ice 225 216 !------------------------------------------------------------------------------! 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) ) ) 217 IF( kt == nit000 ) THEN 218 ! ocean/land mask 219 ALLOCATE( fimask(jpi,jpj) ) 220 IF( rn_ishlat == 0._wp ) THEN 221 DO_2D( 0, 0, 0, 0 ) 222 fimask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) 223 END_2D 224 ELSE 225 DO_2D( 0, 0, 0, 0 ) 226 fimask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) 227 ! Lateral boundary conditions on velocity (modify fimask) 228 IF( fimask(ji,jj) == 0._wp ) THEN 229 fimask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), & 230 & vmask(ji,jj,1), vmask(ji+1,jj,1) ) ) 231 ENDIF 232 END_2D 237 233 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 ) 234 CALL lbc_lnk( 'icedyn_rhg_eap', fimask, 'F', 1.0_wp ) 235 ENDIF 256 236 257 237 !------------------------------------------------------------------------------! … … 401 381 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 382 & + ( 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)383 & ) * r1_e1e2f(ji,jj) * fimask(ji,jj) 404 384 405 385 END_2D … … 760 740 761 741 ! convergence test 762 IF( nn_rhg_chkcvg == 2 ) CALL rhg_cvg_eap( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice )742 IF( nn_rhg_chkcvg == 2 ) CALL rhg_cvg_eap( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice, zmsk15 ) 763 743 ! 764 744 ! ! ==================== ! … … 777 757 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 758 & + ( 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)759 & ) * r1_e1e2f(ji,jj) * fimask(ji,jj) 780 760 781 761 END_2D … … 830 810 & ztauy_ai, 'V', -1.0_wp, ztaux_bi, 'U', -1.0_wp, ztauy_bi, 'V', -1.0_wp ) 831 811 ! 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 )812 CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 ) 813 CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 ) 814 CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 ) 815 CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 ) 816 CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 ) 817 CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 ) 838 818 ENDIF 839 819 840 820 ! --- 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 ) ! strength821 IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * zmsk00 ) ! divergence 822 IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * zmsk00 ) ! shear 823 IF( iom_use('icedlt') ) CALL iom_put( 'icedlt' , pdelta_i * zmsk00 ) ! delta 824 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength 845 825 846 826 ! --- Stress tensor invariants (SIMIP diags) --- ! … … 867 847 ! 868 848 ! 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 stress849 IF( iom_use('normstr') ) CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress 850 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress 871 851 872 852 DEALLOCATE ( zsig_I, zsig_II ) … … 914 894 CALL lbc_lnk( 'icedyn_rhg_eap', zyield11, 'T', 1.0_wp, zyield22, 'T', 1.0_wp, zyield12, 'T', 1.0_wp ) 915 895 916 CALL iom_put( 'yield11', zyield11 * aimsk00 )917 CALL iom_put( 'yield22', zyield22 * aimsk00 )918 CALL iom_put( 'yield12', zyield12 * aimsk00 )896 CALL iom_put( 'yield11', zyield11 * zmsk00 ) 897 CALL iom_put( 'yield22', zyield22 * zmsk00 ) 898 CALL iom_put( 'yield12', zyield12 * zmsk00 ) 919 899 ENDIF 920 900 … … 922 902 IF( iom_use('aniso') ) THEN 923 903 CALL lbc_lnk( 'icedyn_rhg_eap', paniso_11, 'T', 1.0_wp ) 924 CALL iom_put( 'aniso' , paniso_11 * aimsk00 )904 CALL iom_put( 'aniso' , paniso_11 * zmsk00 ) 925 905 ENDIF 926 906 … … 933 913 & zfU, 'U', -1.0_wp, zfV, 'V', -1.0_wp ) 934 914 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)915 CALL iom_put( 'dssh_dx' , zspgU * zmsk00 ) ! Sea-surface tilt term in force balance (x) 916 CALL iom_put( 'dssh_dy' , zspgV * zmsk00 ) ! Sea-surface tilt term in force balance (y) 917 CALL iom_put( 'corstrx' , zCorU * zmsk00 ) ! Coriolis force term in force balance (x) 918 CALL iom_put( 'corstry' , zCorV * zmsk00 ) ! Coriolis force term in force balance (y) 919 CALL iom_put( 'intstrx' , zfU * zmsk00 ) ! Internal force term in force balance (x) 920 CALL iom_put( 'intstry' , zfV * zmsk00 ) ! Internal force term in force balance (y) 941 921 ENDIF 942 922 … … 949 929 DO_2D( 0, 0, 0, 0 ) 950 930 ! 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)931 zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj) 932 zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj) 953 933 954 934 zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component … … 984 964 IF( ln_aEVP ) THEN ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 985 965 CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , & 986 & ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * aimsk15(:,:) )966 & ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) ) 987 967 ELSE ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 988 968 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(:,:) )969 & ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) ) 990 970 ENDIF 991 971 ENDIF … … 995 975 996 976 997 SUBROUTINE rhg_cvg_eap( kt, kiter, kitermax, pu, pv, pub, pvb )977 SUBROUTINE rhg_cvg_eap( kt, kiter, kitermax, pu, pv, pub, pvb, pmsk15 ) 998 978 !!---------------------------------------------------------------------- 999 979 !! *** ROUTINE rhg_cvg_eap *** … … 1010 990 INTEGER , INTENT(in) :: kt, kiter, kitermax ! ocean time-step index 1011 991 REAL(wp), DIMENSION(:,:), INTENT(in) :: pu, pv, pub, pvb ! now and before velocities 992 REAL(wp), DIMENSION(:,:), INTENT(in) :: pmsk15 1012 993 !! 1013 994 INTEGER :: it, idtime, istatus … … 1044 1025 zresm = 0._wp 1045 1026 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) 1027 zresm = 0._wp 1028 DO_2D( 0, 0, 0, 0 ) 1029 zresm = MAX( zresm, MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 1030 & ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * pmsk15(ji,jj) ) 1049 1031 END_2D 1050 1051 zresm = MAXVAL( eap_res )1052 1032 CALL mpp_max( 'icedyn_rhg_evp', zresm ) ! max over the global domain 1053 1033 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
