- Timestamp:
- 2016-06-27T19:20:57+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
r6584 r6746 9 9 !! 3.3 ! 2009-05 (G.Garric) addition of the lim2_evp cas 10 10 !! 3.4 ! 2011-01 (A. Porter) dynamical allocation 11 !! 3.5 ! 2012-08 (R. Benshila) AGRIF 12 !! ----------------------------------------------------------------------13 #if defined key_lim3 || ( defined key_lim2 && ! defined key_lim2_vp ) 14 !!---------------------------------------------------------------------- 15 !! 'key_lim3' OR LIM-3 sea-ice model16 !! 'key_lim 2' AND NOT 'key_lim2_vp' EVP LIM-2sea-ice model11 !! 3.5 ! 2012-08 (R. Benshila) AGRIF 12 !! 3.6 ! 2016-06 (C. Rousset) Rewriting and add the possibility to use mEVP from Bouillon 2013 13 !!---------------------------------------------------------------------- 14 #if defined key_lim3 15 !!---------------------------------------------------------------------- 16 !! 'key_lim3' LIM-3 sea-ice model 17 17 !!---------------------------------------------------------------------- 18 18 !! lim_rhg : computes ice velocities … … 24 24 USE sbc_oce ! Surface boundary condition: ocean fields 25 25 USE sbc_ice ! Surface boundary condition: ice fields 26 #if defined key_lim3 27 USE ice ! LIM-3: ice variables 28 USE dom_ice ! LIM-3: ice domain 29 USE limitd_me ! LIM-3: 30 #else 31 USE ice_2 ! LIM-2: ice variables 32 USE dom_ice_2 ! LIM-2: ice domain 33 #endif 26 USE ice ! ice variables 27 USE limitd_me ! ice strength 34 28 USE lbclnk ! Lateral Boundary Condition / MPP link 35 29 USE lib_mpp ! MPP library … … 38 32 USE prtctl ! Print control 39 33 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 40 #if defined key_agrif && defined key_lim2 41 USE agrif_lim2_interp 42 #endif 43 #if defined key_agrif && defined key_lim3 34 #if defined key_agrif 44 35 USE agrif_lim3_interp 45 36 #endif … … 51 42 PRIVATE 52 43 53 PUBLIC lim_rhg ! routine called by lim_dyn (or lim_dyn_2)44 PUBLIC lim_rhg ! routine called by lim_dyn 54 45 55 46 !! * Substitutions … … 62 53 CONTAINS 63 54 64 SUBROUTINE lim_rhg ( k_j1, k_jpj )55 SUBROUTINE lim_rhg 65 56 !!------------------------------------------------------------------- 66 57 !! *** SUBROUTINE lim_rhg *** … … 109 100 !! e.g. in the Canadian Archipelago 110 101 !! 102 !! ** Notes : There is the possibility to use mEVP from Bouillon 2013 103 !! (by uncommenting some lines in part 3 and changing alpha and beta parameters) 104 !! but this solution appears very unstable (see Kimmritz et al 2016) 105 !! 111 106 !! References : Hunke and Dukowicz, JPO97 112 107 !! Bouillon et al., Ocean Modelling 2009 108 !! Bouillon et al., Ocean Modelling 2013 113 109 !!------------------------------------------------------------------- 114 INTEGER, INTENT(in) :: k_j1 ! southern j-index for ice computation 115 INTEGER, INTENT(in) :: k_jpj ! northern j-index for ice computation 116 !! 117 INTEGER :: ji, jj ! dummy loop indices 118 INTEGER :: jter ! local integers 110 INTEGER :: ji, jj ! dummy loop indices 111 INTEGER :: jter ! local integers 119 112 CHARACTER (len=50) :: charout 120 REAL(wp) :: zt11, zt12, zt21, zt22, ztagnx, ztagny, delta ! 121 REAL(wp) :: za, zstms ! local scalars 122 REAL(wp) :: zc1, zc2, zc3 ! ice mass 123 124 REAL(wp) :: dtevp , z1_dtevp ! time step for subcycling 125 REAL(wp) :: dtotel, z1_dtotel, ecc2, ecci ! square of yield ellipse eccenticity 126 REAL(wp) :: z0, zr, zcca, zccb ! temporary scalars 127 REAL(wp) :: zu_ice2, zv_ice1 ! 128 REAL(wp) :: zddc, zdtc ! delta on corners and on centre 129 REAL(wp) :: zdst ! shear at the center of the grid point 130 REAL(wp) :: zdsshx, zdsshy ! term for the gradient of ocean surface 131 REAL(wp) :: sigma1, sigma2 ! internal ice stress 132 133 REAL(wp) :: zresm ! Maximal error on ice velocity 134 REAL(wp) :: zintb, zintn ! dummy argument 135 136 REAL(wp), POINTER, DIMENSION(:,:) :: zpresh ! temporary array for ice strength 137 REAL(wp), POINTER, DIMENSION(:,:) :: zpreshc ! Ice strength on grid cell corners (zpreshc) 138 REAL(wp), POINTER, DIMENSION(:,:) :: zfrld1, zfrld2 ! lead fraction on U/V points 139 REAL(wp), POINTER, DIMENSION(:,:) :: zmass1, zmass2 ! ice/snow mass on U/V points 140 REAL(wp), POINTER, DIMENSION(:,:) :: zcorl1, zcorl2 ! coriolis parameter on U/V points 141 REAL(wp), POINTER, DIMENSION(:,:) :: za1ct , za2ct ! temporary arrays 142 REAL(wp), POINTER, DIMENSION(:,:) :: v_oce1 ! ocean u/v component on U points 143 REAL(wp), POINTER, DIMENSION(:,:) :: u_oce2 ! ocean u/v component on V points 144 REAL(wp), POINTER, DIMENSION(:,:) :: u_ice2, v_ice1 ! ice u/v component on V/U point 145 REAL(wp), POINTER, DIMENSION(:,:) :: zf1 , zf2 ! arrays for internal stresses 146 REAL(wp), POINTER, DIMENSION(:,:) :: zmask ! mask ocean grid points 113 114 REAL(wp) :: dtevp, z1_dtevp ! time step for subcycling 115 REAL(wp) :: ecc2, z1_ecc2 ! square of yield ellipse eccenticity 116 REAL(wp) :: zbeta, zalph1, z1_alph1, zalph2, z1_alph2 ! alpha and beta from Bouillon 2009 and 2013 117 REAL(wp) :: zm1, zm2, zm3 ! ice/snow mass 118 REAL(wp) :: delta, zp_delf, zds2 119 REAL(wp) :: zTauO, zTauA, zTauB, ZCor, zmt, zu_ice2, zv_ice1, zvel ! temporary scalars 120 121 REAL(wp) :: zsig1, zsig2 ! internal ice stress 122 REAL(wp) :: zresm ! Maximal error on ice velocity 123 REAL(wp) :: zintb, zintn ! dummy argument 124 125 REAL(wp), POINTER, DIMENSION(:,:) :: zpresh ! temporary array for ice strength 126 REAL(wp), POINTER, DIMENSION(:,:) :: z1_e1t0, z1_e2t0 ! scale factors 127 REAL(wp), POINTER, DIMENSION(:,:) :: zp_delt ! P/delta at T points 128 ! 129 REAL(wp), POINTER, DIMENSION(:,:) :: za1 , za2 ! ice fraction on U/V points 130 REAL(wp), POINTER, DIMENSION(:,:) :: zmass1, zmass2 ! ice/snow mass on U/V points 131 REAL(wp), POINTER, DIMENSION(:,:) :: zcor1 , zcor2 ! coriolis parameter on U/V points 132 REAL(wp), POINTER, DIMENSION(:,:) :: zspgu , zspgv ! surface pressure gradient at U/V points 133 REAL(wp), POINTER, DIMENSION(:,:) :: v_oce1, u_oce2, v_ice1, u_ice2 ! ocean/ice u/v component on V/U points 134 REAL(wp), POINTER, DIMENSION(:,:) :: zf1 , zf2 ! internal stresses 147 135 148 REAL(wp), POINTER, DIMENSION(:,:) :: zdt ! tension at centre of grid cells 149 REAL(wp), POINTER, DIMENSION(:,:) :: zds ! Shear on northeast corner of grid cells 150 REAL(wp), POINTER, DIMENSION(:,:) :: zs1 , zs2 ! Diagonal stress tensor components zs1 and zs2 151 REAL(wp), POINTER, DIMENSION(:,:) :: zs12 ! Non-diagonal stress tensor component zs12 152 REAL(wp), POINTER, DIMENSION(:,:) :: zu_ice, zv_ice, zresr ! Local error on velocity 153 REAL(wp), POINTER, DIMENSION(:,:) :: zpice ! array used for the calculation of ice surface slope: 154 ! ocean surface (ssh_m) if ice is not embedded 155 ! ice top surface if ice is embedded 156 157 REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter 158 REAL(wp), PARAMETER :: zvmin = 1.0e-03_wp ! ice volume below which ice velocity equals ocean velocity 136 REAL(wp), POINTER, DIMENSION(:,:) :: zdt, zds ! tension and shear 137 REAL(wp), POINTER, DIMENSION(:,:) :: zs1, zs2, zs12 ! stress tensor components 138 REAL(wp), POINTER, DIMENSION(:,:) :: zu_ice, zv_ice, zresr ! check convergence 139 REAL(wp), POINTER, DIMENSION(:,:) :: zpice ! array used for the calculation of ice surface slope: 140 ! ocean surface (ssh_m) if ice is not embedded 141 ! ice top surface if ice is embedded 142 REAL(wp), POINTER, DIMENSION(:,:) :: zswitch1, zswitch2 ! dummy arrays 143 144 REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter 145 REAL(wp), PARAMETER :: zvmin = 1.0e-03_wp ! ice volume below which ice velocity equals ocean velocity 159 146 !!------------------------------------------------------------------- 160 147 161 CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 162 CALL wrk_alloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask ) 163 CALL wrk_alloc( jpi,jpj, zf1 , zu_ice, zf2 , zv_ice , zdt , zds ) 164 CALL wrk_alloc( jpi,jpj, zs1 , zs2 , zs12 , zresr , zpice ) 165 166 #if defined key_lim2 && ! defined key_lim2_vp 167 # if defined key_agrif 168 USE ice_2, vt_s => hsnm 169 USE ice_2, vt_i => hicm 170 # else 171 vt_s => hsnm 172 vt_i => hicm 173 # endif 174 at_i(:,:) = 1. - frld(:,:) 175 #endif 176 #if defined key_agrif && defined key_lim2 177 CALL agrif_rhg_lim2_load ! First interpolation of coarse values 178 #endif 179 #if defined key_agrif && defined key_lim3 148 CALL wrk_alloc( jpi,jpj, zpresh, z1_e1t0, z1_e2t0, zp_delt ) 149 CALL wrk_alloc( jpi,jpj, za1, za2, zmass1, zmass2, zcor1, zcor2 ) 150 CALL wrk_alloc( jpi,jpj, zspgu, zspgv, v_oce1, u_oce2, v_ice1, u_ice2, zf1, zf2 ) 151 CALL wrk_alloc( jpi,jpj, zds, zdt, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 152 CALL wrk_alloc( jpi,jpj, zswitch1, zswitch2 ) 153 154 #if defined key_agrif 180 155 CALL agrif_interp_lim3('U') ! First interpolation of coarse values 181 CALL agrif_interp_lim3('V') ! First interpolation of coarse values 182 #endif 183 ! 184 !------------------------------------------------------------------------------! 185 ! 1) Ice strength (zpresh) ! 186 !------------------------------------------------------------------------------! 187 ! 188 ! Put every vector to 0 189 delta_i(:,:) = 0._wp ; 190 zpresh (:,:) = 0._wp ; 191 zpreshc(:,:) = 0._wp 192 u_ice2 (:,:) = 0._wp ; v_ice1(:,:) = 0._wp 193 divu_i (:,:) = 0._wp ; zdt (:,:) = 0._wp ; zds(:,:) = 0._wp 194 shear_i(:,:) = 0._wp 195 196 #if defined key_lim3 197 CALL lim_itd_me_icestrength( nn_icestr ) ! LIM-3: Ice strength on T-points 198 #endif 199 200 DO jj = k_j1 , k_jpj ! Ice mass and temp variables 201 DO ji = 1 , jpi 202 #if defined key_lim3 203 zpresh(ji,jj) = tmask(ji,jj,1) * strength(ji,jj) 204 #endif 205 #if defined key_lim2 206 zpresh(ji,jj) = tmask(ji,jj,1) * pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) ) 207 #endif 208 ! zmask = 1 where there is ice or on land 209 zmask(ji,jj) = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - zepsi ) ) ) * tmask(ji,jj,1) 156 CALL agrif_interp_lim3('V') 157 #endif 158 ! 159 !------------------------------------------------------------------------------! 160 ! 0) define some variables 161 !------------------------------------------------------------------------------! 162 ! ecc2: square of yield ellipse eccenticrity 163 ecc2 = rn_ecc * rn_ecc 164 z1_ecc2 = 1._wp / ecc2 165 166 ! Time step for subcycling 167 dtevp = rdt_ice / REAL( nn_nevp ) 168 z1_dtevp = 1._wp / dtevp 169 170 ! alpha parameters (Bouillon 2009) 171 zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp 172 zalph2 = zalph1 * z1_ecc2 173 174 ! alpha and beta parameters (Bouillon 2013) 175 !!zalph1 = 40. 176 !!zalph2 = 40. 177 !!zbeta = 3000. 178 !!zbeta = REAL( nn_nevp ) ! close to classical EVP of Hunke (2001) 179 180 z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 181 z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 182 183 !------------------------------------------------------------------------------! 184 ! 1) initialize arrays 185 !------------------------------------------------------------------------------! 186 ! Initialise stress tensor 187 zs1 (:,:) = stress1_i (:,:) 188 zs2 (:,:) = stress2_i (:,:) 189 zs12(:,:) = stress12_i(:,:) 190 191 ! Ice strength 192 CALL lim_itd_me_icestrength( nn_icestr ) 193 zpresh(:,:) = tmask(:,:,1) * strength(:,:) 194 195 ! scale factors 196 DO jj = 2, jpjm1 197 DO ji = fs_2, fs_jpim1 198 z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj ) + e1t(ji,jj ) ) 199 z1_e2t0(ji,jj) = 1._wp / ( e2t(ji ,jj+1) + e2t(ji,jj ) ) 210 200 END DO 211 201 END DO 212 213 ! Ice strength on grid cell corners (zpreshc) 214 ! needed for calculation of shear stress 215 DO jj = k_j1+1, k_jpj-1 216 DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1) 217 zstms = tmask(ji+1,jj+1,1) * wght(ji+1,jj+1,2,2) + tmask(ji,jj+1,1) * wght(ji+1,jj+1,1,2) + & 218 & tmask(ji+1,jj,1) * wght(ji+1,jj+1,2,1) + tmask(ji,jj,1) * wght(ji+1,jj+1,1,1) 219 zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) + & 220 & zpresh(ji+1,jj) * wght(ji+1,jj+1,2,1) + zpresh(ji,jj) * wght(ji+1,jj+1,1,1) & 221 & ) / MAX( zstms, zepsi ) 222 END DO 223 END DO 224 CALL lbc_lnk( zpreshc(:,:), 'F', 1. ) 202 225 203 ! 226 204 !------------------------------------------------------------------------------! 227 205 ! 2) Wind / ocean stress, mass terms, coriolis terms 228 206 !------------------------------------------------------------------------------! 229 !230 ! Wind stress, coriolis and mass terms on the sides of the squares231 ! zfrld1: lead fraction on U-points232 ! zfrld2: lead fraction on V-points233 ! zmass1: ice/snow mass on U-points234 ! zmass2: ice/snow mass on V-points235 ! zcorl1: Coriolis parameter on U-points236 ! zcorl2: Coriolis parameter on V-points237 ! (ztagnx,ztagny): wind stress on U/V points238 ! v_oce1: ocean v component on u points239 ! u_oce2: ocean u component on v points240 207 241 208 IF( nn_ice_embd == 2 ) THEN !== embedded sea ice: compute representative ice top surface ==! … … 249 216 zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 250 217 ! 251 zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:)) * r1_rau0218 zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0 252 219 ! 253 220 ELSE !== non-embedded sea ice: use ocean surface for slope calculation ==! … … 255 222 ENDIF 256 223 257 DO jj = k_j1+1, k_jpj-1224 DO jj = 2, jpjm1 258 225 DO ji = fs_2, fs_jpim1 259 226 260 zc1 = tmask(ji ,jj ,1) * ( rhosn * vt_s(ji ,jj ) + rhoic * vt_i(ji ,jj ) ) 261 zc2 = tmask(ji+1,jj ,1) * ( rhosn * vt_s(ji+1,jj ) + rhoic * vt_i(ji+1,jj ) ) 262 zc3 = tmask(ji ,jj+1,1) * ( rhosn * vt_s(ji ,jj+1) + rhoic * vt_i(ji ,jj+1) ) 263 264 zt11 = tmask(ji ,jj,1) * e1t(ji ,jj) 265 zt12 = tmask(ji+1,jj,1) * e1t(ji+1,jj) 266 zt21 = tmask(ji,jj ,1) * e2t(ji,jj ) 267 zt22 = tmask(ji,jj+1,1) * e2t(ji,jj+1) 268 269 ! Leads area. 270 zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + zepsi ) 271 zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + zepsi ) 272 273 ! Mass, coriolis coeff. and currents 274 zmass1(ji,jj) = ( zt12 * zc1 + zt11 * zc2 ) / ( zt11 + zt12 + zepsi ) 275 zmass2(ji,jj) = ( zt22 * zc1 + zt21 * zc3 ) / ( zt21 + zt22 + zepsi ) 276 zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj) * fcor(ji,jj) + e1t(ji,jj) * fcor(ji+1,jj) ) & 277 & / ( e1t(ji,jj) + e1t(ji+1,jj) + zepsi ) 278 zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1) * fcor(ji,jj) + e2t(ji,jj) * fcor(ji,jj+1) ) & 279 & / ( e2t(ji,jj+1) + e2t(ji,jj) + zepsi ) 280 ! 281 ! Ocean has no slip boundary condition 282 v_oce1(ji,jj) = 0.5 * ( ( v_oce(ji ,jj) + v_oce(ji ,jj-1) ) * e1t(ji,jj) & 283 & + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji+1,jj) ) & 284 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 285 286 u_oce2(ji,jj) = 0.5 * ( ( u_oce(ji,jj ) + u_oce(ji-1,jj ) ) * e2t(ji,jj) & 287 & + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj+1) ) & 288 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 289 290 ! Wind stress at U,V-point 291 ztagnx = ( 1. - zfrld1(ji,jj) ) * utau_ice(ji,jj) 292 ztagny = ( 1. - zfrld2(ji,jj) ) * vtau_ice(ji,jj) 293 294 ! Computation of the velocity field taking into account the ice internal interaction. 295 ! Terms that are independent of the velocity field. 296 297 ! SB On utilise maintenant le gradient de la pente de l'ocean 298 ! include it later 299 300 zdsshx = ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj) 301 zdsshy = ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj) 302 303 za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx 304 za2ct(ji,jj) = ztagny - zmass2(ji,jj) * grav * zdsshy 227 ! ice fraction at U-V points 228 za1(ji,jj) = ( at_i(ji,jj) * e1t(ji+1,jj) + at_i(ji+1,jj) * e1t(ji,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 229 za2(ji,jj) = ( at_i(ji,jj) * e2t(ji,jj+1) + at_i(ji,jj+1) * e2t(ji,jj) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 230 231 ! Ice/snow mass at U-V points 232 zm1 = ( rhosn * vt_s(ji ,jj ) + rhoic * vt_i(ji ,jj ) ) 233 zm2 = ( rhosn * vt_s(ji+1,jj ) + rhoic * vt_i(ji+1,jj ) ) 234 zm3 = ( rhosn * vt_s(ji ,jj+1) + rhoic * vt_i(ji ,jj+1) ) 235 zmass1(ji,jj) = ( zm1 * e1t(ji+1,jj) + zm2 * e1t(ji,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 236 zmass2(ji,jj) = ( zm1 * e2t(ji,jj+1) + zm3 * e2t(ji,jj) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 237 238 ! Ocean currents at U-V points 239 v_oce1(ji,jj) = 0.5_wp * ( ( v_oce(ji ,jj) + v_oce(ji ,jj-1) ) * e1t(ji+1,jj) & 240 & + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 241 242 u_oce2(ji,jj) = 0.5_wp * ( ( u_oce(ji,jj ) + u_oce(ji-1,jj ) ) * e2t(ji,jj+1) & 243 & + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 244 245 ! Coriolis at U-V points (m*f) 246 zcor1(ji,jj) = zmass1(ji,jj) * 0.5_wp * ( ff(ji,jj) + ff(ji ,jj-1) ) 247 zcor2(ji,jj) = - zmass2(ji,jj) * 0.5_wp * ( ff(ji,jj) + ff(ji-1,jj ) ) 248 249 ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points 250 zspgu(ji,jj) = - zmass1(ji,jj) * grav * ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj) 251 zspgv(ji,jj) = - zmass2(ji,jj) * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj) 252 253 ! switches 254 zswitch1(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) 255 zswitch2(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) 305 256 306 257 END DO … … 312 263 !------------------------------------------------------------------------------! 313 264 ! 314 ! Time step for subcycling315 dtevp = rdt_ice / nn_nevp316 #if defined key_lim3317 dtotel = dtevp / ( 2._wp * rn_relast * rdt_ice )318 #else319 dtotel = dtevp / ( 2._wp * telast )320 #endif321 z1_dtotel = 1._wp / ( 1._wp + dtotel )322 z1_dtevp = 1._wp / dtevp323 !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter)324 ecc2 = rn_ecc * rn_ecc325 ecci = 1. / ecc2326 327 !-Initialise stress tensor328 zs1 (:,:) = stress1_i (:,:)329 zs2 (:,:) = stress2_i (:,:)330 zs12(:,:) = stress12_i(:,:)331 332 265 ! !----------------------! 333 266 DO jter = 1 , nn_nevp ! loop over jter ! 334 267 ! !----------------------! 335 DO jj = k_j1, k_jpj-1 336 zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 337 zv_ice(:,jj) = v_ice(:,jj) 338 END DO 339 340 DO jj = k_j1+1, k_jpj-1 341 DO ji = fs_2, fs_jpim1 !RB bug no vect opt due to zmask 342 343 ! 344 !- Divergence, tension and shear (Section a. Appendix B of Hunke & Dukowicz, 2002) 345 !- divu_i(:,:), zdt(:,:): divergence and tension at centre of grid cells 346 !- zds(:,:): shear on northeast corner of grid cells 347 ! 348 !- IMPORTANT REMINDER: Dear Gurvan, note that, the way these terms are coded, 349 ! there are many repeated calculations. 350 ! Speed could be improved by regrouping terms. For 351 ! the moment, however, the stress is on clarity of coding to avoid 352 ! bugs (Martin, for Miguel). 353 ! 354 !- ALSO: arrays zdt, zds and delta could 355 ! be removed in the future to minimise memory demand. 356 ! 357 !- MORE NOTES: Note that we are calculating deformation rates and stresses on the corners of 358 ! grid cells, exactly as in the B grid case. For simplicity, the indexation on 359 ! the corners is the same as in the B grid. 360 ! 361 ! 362 divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 363 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 268 ! Convergence test 269 IF(ln_ctl) THEN 270 DO jj = 1, jpjm1 271 zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 272 zv_ice(:,jj) = v_ice(:,jj) 273 END DO 274 ENDIF 275 276 ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 277 DO jj = 2, jpjm1 278 DO ji = fs_2, fs_jpim1 279 280 ! divergence at T points 281 divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 282 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 364 283 & ) * r1_e12t(ji,jj) 365 284 285 ! tension at T points 366 286 zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & 367 287 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 368 288 & ) * r1_e12t(ji,jj) 369 289 370 ! 290 ! shear at F points 371 291 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) & 372 292 & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 373 & ) * r1_e12f(ji,jj) * ( 2._wp - fmask(ji,jj,1) ) & 374 & * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1) 375 376 377 v_ice1(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji+1,jj) & 378 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) & 379 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 380 381 u_ice2(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) & 382 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) & 383 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 384 END DO 385 END DO 386 387 CALL lbc_lnk_multi( v_ice1, 'U', -1., u_ice2, 'V', -1. ) ! lateral boundary cond. 293 & ) * r1_e12f(ji,jj) 294 295 ! u_ice at V point 296 u_ice2(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) & 297 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 298 299 ! v_ice at U point 300 v_ice1(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji+1,jj) & 301 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 302 303 END DO 304 END DO 305 CALL lbc_lnk( zds, 'F', 1. ) 306 307 DO jj = 2, jpjm1 308 DO ji = 2, jpim1 ! no vector loop 309 310 ! shear**2 at T points (doc eq. A16) 311 zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e12f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e12f(ji-1,jj ) & 312 & + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1) & 313 & ) * 0.25_wp * r1_e12t(ji,jj) 314 315 ! delta at T points 316 delta = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt(ji,jj) * zdt(ji,jj) + zds2 ) * usecc2 ) 317 delta_i(ji,jj) = delta + rn_creepl 318 319 ! P/delta at T points 320 zp_delt(ji,jj) = zpresh(ji,jj) / delta_i(ji,jj) 321 END DO 322 END DO 323 CALL lbc_lnk( zp_delt, 'T', 1. ) 324 325 ! --- Stress tensor --- ! 326 DO jj = 2, jpjm1 327 DO ji = 2, jpim1 ! no vector loop 328 329 ! P/delta at F points 330 zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) ) 331 332 ! stress tensor at T points 333 zs1(ji,jj) = ( zs1 (ji,jj) * zalph1 + zp_delt(ji,jj) * ( divu_i(ji,jj) - (delta_i(ji,jj)-rn_creepl) ) ) * z1_alph1 334 zs2(ji,jj) = ( zs2 (ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt(ji,jj) * z1_ecc2 ) ) * z1_alph2 335 ! F points 336 zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2 337 END DO 338 END DO 339 CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 340 341 ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- ! 342 DO jj = 2, jpjm1 343 DO ji = fs_2, fs_jpim1 344 ! U points 345 zf1(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) & 346 & + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj) & 347 & ) * r1_e2u(ji,jj) & 348 & + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) & 349 & ) * 2._wp * r1_e1u(ji,jj) & 350 & ) * r1_e12u(ji,jj) 351 352 ! V points 353 zf2(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) & 354 & - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) & 355 & ) * r1_e1v(ji,jj) & 356 & + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) & 357 & ) * 2._wp * r1_e2v(ji,jj) & 358 & ) * r1_e12v(ji,jj) 359 END DO 360 END DO 361 ! 362 ! --- Computation of ice velocity --- ! 363 ! Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen smart and large nn_nevp 364 ! Bouillon et al. 2009 (eq 34-35) => stable 365 IF( MOD(jter,2) .EQ. 0 ) THEN ! even iterations 366 367 DO jj = 2, jpjm1 368 DO ji = fs_2, fs_jpim1 369 370 zvel = SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) & 371 & + ( u_ice2(ji,jj) - u_oce2(ji,jj) ) * ( u_ice2(ji,jj) - u_oce2(ji,jj) ) ) 372 373 zTauA = za2(ji,jj) * vtau_ice(ji,jj) 374 zTauO = za2(ji,jj) * rhoco * zvel 375 ztauB = - tau_icebfr(ji,jj) / zvel 376 zCor = zcor2(ji,jj) * u_ice2(ji,jj) 377 zmt = zmass2(ji,jj) * z1_dtevp 378 379 ! Bouillon 2009 380 v_ice(ji,jj) = ( zmt * v_ice(ji,jj) + zf2(ji,jj) + zCor + zTauA + zTauO * v_oce(ji,jj) + zspgv(ji,jj) & 381 & ) / MAX( zmt + zTauO - zTauB, zepsi ) * zswitch2(ji,jj) 382 ! Bouillon 2013 383 !!v_ice(ji,jj) = ( zmt * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) ) & 384 !! & + zf2(ji,jj) + zCor + zTauA + zTauO * v_oce(ji,jj) + zspgv(ji,jj) & 385 !! & ) / MAX( zmt * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch2(ji,jj) 386 387 END DO 388 END DO 389 CALL lbc_lnk( v_ice, 'V', -1. ) 390 391 #if defined key_agrif 392 CALL agrif_interp_lim3( 'V' ) 393 #endif 394 #if defined key_bdy 395 CALL bdy_ice_lim_dyn( 'V' ) 396 #endif 397 398 DO jj = 2, jpjm1 399 DO ji = fs_2, fs_jpim1 400 401 ! v_ice at U point 402 zv_ice1 = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji+1,jj) & 403 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 404 405 zvel = SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) & 406 & + ( v_ice1(ji,jj) - v_oce1(ji,jj) ) * ( v_ice1(ji,jj) - v_oce1(ji,jj) ) ) 407 408 zTauA = za1(ji,jj) * utau_ice(ji,jj) 409 zTauO = za1(ji,jj) * rhoco * zvel 410 ztauB = - tau_icebfr(ji,jj) / zvel 411 zCor = zcor1(ji,jj) * zv_ice1 412 zmt = zmass1(ji,jj) * z1_dtevp 413 414 ! Bouillon 2009 415 u_ice(ji,jj) = ( zmt * u_ice(ji,jj) + zf1(ji,jj) + zCor + zTauA + zTauO * u_oce(ji,jj) + zspgu(ji,jj) & 416 & ) / MAX( zmt + zTauO - zTauB, zepsi ) * zswitch1(ji,jj) 417 ! Bouillon 2013 418 !!u_ice(ji,jj) = ( zmt * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) ) & 419 !! & + zf1(ji,jj) + zCor + zTauA + zTauO * u_oce(ji,jj) + zspgu(ji,jj) & 420 !! & ) / MAX( zmt * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch1(ji,jj) 421 END DO 422 END DO 423 CALL lbc_lnk( u_ice, 'U', -1. ) 424 425 #if defined key_agrif 426 CALL agrif_interp_lim3( 'U' ) 427 #endif 428 #if defined key_bdy 429 CALL bdy_ice_lim_dyn( 'U' ) 430 #endif 431 432 ELSE ! odd iterations 433 434 DO jj = 2, jpjm1 435 DO ji = fs_2, fs_jpim1 436 437 zvel = SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) & 438 & + ( v_ice1(ji,jj) - v_oce1(ji,jj) ) * ( v_ice1(ji,jj) - v_oce1(ji,jj) ) ) 439 440 zTauA = za1(ji,jj) * utau_ice(ji,jj) 441 zTauO = za1(ji,jj) * rhoco * zvel 442 ztauB = - tau_icebfr(ji,jj) / zvel 443 zCor = zcor1(ji,jj) * v_ice1(ji,jj) 444 zmt = zmass1(ji,jj) * z1_dtevp 445 446 ! Bouillon 2009 447 u_ice(ji,jj) = ( zmt * u_ice(ji,jj) + zf1(ji,jj) + zCor + zTauA + zTauO * u_oce(ji,jj) + zspgu(ji,jj) & 448 & ) / MAX( zmt + zTauO - zTauB, zepsi ) * zswitch1(ji,jj) 449 ! Bouillon 2013 450 !!u_ice(ji,jj) = ( zmt * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) ) & 451 !! & + zf1(ji,jj) + zCor + zTauA + zTauO * u_oce(ji,jj) + zspgu(ji,jj) & 452 !! & ) / MAX( zmt * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch1(ji,jj) 453 END DO 454 END DO 455 CALL lbc_lnk( u_ice, 'U', -1. ) 456 457 #if defined key_agrif 458 CALL agrif_interp_lim3( 'U' ) 459 #endif 460 #if defined key_bdy 461 CALL bdy_ice_lim_dyn( 'U' ) 462 #endif 463 464 DO jj = 2, jpjm1 465 DO ji = fs_2, fs_jpim1 466 467 ! u_ice at V point 468 zu_ice2 = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) & 469 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 470 471 zvel = SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) & 472 & + ( u_ice2(ji,jj) - u_oce2(ji,jj) ) * ( u_ice2(ji,jj) - u_oce2(ji,jj) ) ) 388 473 389 DO jj = k_j1+1, k_jpj-1 390 DO ji = fs_2, fs_jpim1 391 392 !- Calculate Delta at centre of grid cells 393 zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj ) * v_ice1(ji-1,jj ) & 394 & + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji ,jj-1) * u_ice2(ji ,jj-1) & 395 & ) * r1_e12t(ji,jj) 396 397 delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 ) 398 delta_i(ji,jj) = delta + rn_creepl 399 400 !- Calculate Delta on corners 401 zddc = ( ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 402 & + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 403 & ) * r1_e12f(ji,jj) 404 405 zdtc = (- ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 406 & + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 407 & ) * r1_e12f(ji,jj) 408 409 zddc = SQRT( zddc**2 + ( zdtc**2 + zds(ji,jj)**2 ) * usecc2 ) + rn_creepl 410 411 !-Calculate stress tensor components zs1 and zs2 at centre of grid cells (see section 3.5 of CICE user's guide). 412 zs1(ji,jj) = ( zs1 (ji,jj) + dtotel * ( divu_i(ji,jj) - delta ) / delta_i(ji,jj) * zpresh(ji,jj) & 413 & ) * z1_dtotel 414 zs2(ji,jj) = ( zs2 (ji,jj) + dtotel * ecci * zdt(ji,jj) / delta_i(ji,jj) * zpresh(ji,jj) & 415 & ) * z1_dtotel 416 !-Calculate stress tensor component zs12 at corners 417 zs12(ji,jj) = ( zs12(ji,jj) + dtotel * ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj) & 418 & ) * z1_dtotel 419 420 END DO 421 END DO 422 423 CALL lbc_lnk_multi( zs1 , 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 424 425 ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) 426 DO jj = k_j1+1, k_jpj-1 427 DO ji = fs_2, fs_jpim1 428 !- contribution of zs1, zs2 and zs12 to zf1 429 zf1(ji,jj) = 0.5 * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) & 430 & + ( zs2(ji+1,jj) * e2t(ji+1,jj)**2 - zs2(ji,jj) * e2t(ji,jj)**2 ) * r1_e2u(ji,jj) & 431 & + 2.0 * ( zs12(ji,jj) * e1f(ji,jj)**2 - zs12(ji,jj-1) * e1f(ji,jj-1)**2 ) * r1_e1u(ji,jj) & 432 & ) * r1_e12u(ji,jj) 433 ! contribution of zs1, zs2 and zs12 to zf2 434 zf2(ji,jj) = 0.5 * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) & 435 & - ( zs2(ji,jj+1) * e1t(ji,jj+1)**2 - zs2(ji,jj) * e1t(ji,jj)**2 ) * r1_e1v(ji,jj) & 436 & + 2.0 * ( zs12(ji,jj) * e2f(ji,jj)**2 - zs12(ji-1,jj) * e2f(ji-1,jj)**2 ) * r1_e2v(ji,jj) & 437 & ) * r1_e12v(ji,jj) 438 END DO 439 END DO 440 ! 441 ! Computation of ice velocity 442 ! 443 ! Both the Coriolis term and the ice-ocean drag are solved semi-implicitly. 444 ! 445 IF (MOD(jter,2).eq.0) THEN 446 447 DO jj = k_j1+1, k_jpj-1 448 DO ji = fs_2, fs_jpim1 449 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1) 450 z0 = zmass1(ji,jj) * z1_dtevp 451 452 ! SB modif because ocean has no slip boundary condition 453 zv_ice1 = 0.5 * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji ,jj) & 454 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) ) & 455 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 456 za = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 + & 457 & ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) ) 458 zr = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj) 459 zcca = z0 + za 460 zccb = zcorl1(ji,jj) 461 u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch 474 zTauA = za2(ji,jj) * vtau_ice(ji,jj) 475 zTauO = za2(ji,jj) * rhoco * zvel 476 ztauB = - tau_icebfr(ji,jj) / zvel 477 zCor = zcor2(ji,jj) * zu_ice2 478 zmt = zmass2(ji,jj) * z1_dtevp 479 480 ! Bouillon 2009 481 v_ice(ji,jj) = ( zmt * v_ice(ji,jj) + zf2(ji,jj) + zCor + zTauA + zTauO * v_oce(ji,jj) + zspgv(ji,jj) & 482 & ) / MAX( zmt + zTauO - zTauB, zepsi ) * zswitch2(ji,jj) 483 ! Bouillon 2013 484 !!v_ice(ji,jj) = ( zmt * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) ) & 485 !! & + zf2(ji,jj) + zCor + zTauA + zTauO * v_oce(ji,jj) + zspgv(ji,jj) & 486 !! & ) / MAX( zmt * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch2(ji,jj) 487 462 488 END DO 463 489 END DO 464 465 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 466 #if defined key_agrif && defined key_lim2 467 CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 468 #endif 469 #if defined key_agrif && defined key_lim3 470 CALL agrif_interp_lim3( 'U' ) 490 CALL lbc_lnk( v_ice, 'V', -1. ) 491 492 #if defined key_agrif 493 CALL agrif_interp_lim3( 'V' ) 471 494 #endif 472 495 #if defined key_bdy 473 CALL bdy_ice_lim_dyn( 'U' ) 474 #endif 475 476 DO jj = k_j1+1, k_jpj-1 477 DO ji = fs_2, fs_jpim1 478 479 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1) 480 z0 = zmass2(ji,jj) * z1_dtevp 481 ! SB modif because ocean has no slip boundary condition 482 zu_ice2 = 0.5 * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj) & 483 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) ) & 484 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 485 za = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 + & 486 & ( v_ice(ji,jj) - v_oce(ji,jj))**2 ) * ( 1.0 - zfrld2(ji,jj) ) 487 zr = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj) 488 zcca = z0 + za 489 zccb = zcorl2(ji,jj) 490 v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch 491 END DO 492 END DO 493 494 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 495 #if defined key_agrif && defined key_lim2 496 CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 497 #endif 498 #if defined key_agrif && defined key_lim3 499 CALL agrif_interp_lim3( 'V' ) 500 #endif 501 #if defined key_bdy 502 CALL bdy_ice_lim_dyn( 'V' ) 503 #endif 504 505 ELSE 506 DO jj = k_j1+1, k_jpj-1 507 DO ji = fs_2, fs_jpim1 508 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1) 509 z0 = zmass2(ji,jj) * z1_dtevp 510 ! SB modif because ocean has no slip boundary condition 511 zu_ice2 = 0.5 * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj) & 512 & +( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) ) & 513 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 514 515 za = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 + & 516 & ( v_ice(ji,jj) - v_oce(ji,jj) )**2 ) * ( 1.0 - zfrld2(ji,jj) ) 517 zr = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj) 518 zcca = z0 + za 519 zccb = zcorl2(ji,jj) 520 v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch 521 END DO 522 END DO 523 524 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 525 #if defined key_agrif && defined key_lim2 526 CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 527 #endif 528 #if defined key_agrif && defined key_lim3 529 CALL agrif_interp_lim3( 'V' ) 530 #endif 531 #if defined key_bdy 532 CALL bdy_ice_lim_dyn( 'V' ) 533 #endif 534 535 DO jj = k_j1+1, k_jpj-1 536 DO ji = fs_2, fs_jpim1 537 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1) 538 z0 = zmass1(ji,jj) * z1_dtevp 539 zv_ice1 = 0.5 * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji,jj) & 540 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) ) & 541 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 542 543 za = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 + & 544 & ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) ) 545 zr = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj) 546 zcca = z0 + za 547 zccb = zcorl1(ji,jj) 548 u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch 549 END DO 550 END DO 551 552 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 553 #if defined key_agrif && defined key_lim2 554 CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 555 #endif 556 #if defined key_agrif && defined key_lim3 557 CALL agrif_interp_lim3( 'U' ) 558 #endif 559 #if defined key_bdy 560 CALL bdy_ice_lim_dyn( 'U' ) 496 CALL bdy_ice_lim_dyn( 'V' ) 561 497 #endif 562 498 563 499 ENDIF 564 500 501 ! Convergence test 565 502 IF(ln_ctl) THEN 566 !--- Convergence test. 567 DO jj = k_j1+1 , k_jpj-1 503 DO jj = 2 , jpjm1 568 504 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 569 505 END DO 570 zresm = MAXVAL( zresr( 1:jpi, k_j1+1:k_jpj-1 ) )506 zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) ) 571 507 IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain 572 508 ENDIF 573 509 ! 574 510 ! ! ==================== ! 575 511 END DO ! end loop over jter ! … … 577 513 ! 578 514 !------------------------------------------------------------------------------! 579 ! 4) Prevent ice velocities when theice is thin515 ! 4) Limit ice velocities when ice is thin 580 516 !------------------------------------------------------------------------------! 581 517 ! If the ice volume is below zvmin then ice velocity should equal the 582 518 ! ocean velocity. This prevents high velocity when ice is thin 583 DO jj = k_j1+1, k_jpj-1519 DO jj = 2, jpjm1 584 520 DO ji = fs_2, fs_jpim1 585 IF ( vt_i(ji,jj) <= zvmin ) THEN 586 u_ice(ji,jj) = u_oce(ji,jj) 587 v_ice(ji,jj) = v_oce(ji,jj) 588 ENDIF 521 rswitch = MAX( 0._wp, SIGN( 1._wp, vt_i(ji,jj) - zvmin ) ) 522 u_ice(ji,jj) = ( rswitch * u_ice(ji,jj) + (1._wp - rswitch) * u_oce(ji,jj) ) * zswitch1(ji,jj) 523 v_ice(ji,jj) = ( rswitch * v_ice(ji,jj) + (1._wp - rswitch) * v_oce(ji,jj) ) * zswitch2(ji,jj) 589 524 END DO 590 525 END DO 591 526 592 CALL lbc_lnk_multi( u_ice(:,:), 'U', -1., v_ice(:,:), 'V', -1. ) 593 594 #if defined key_agrif && defined key_lim2 595 CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'U' ) 596 CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'V' ) 597 #endif 598 #if defined key_agrif && defined key_lim3 527 CALL lbc_lnk_multi( u_ice, 'U', -1., v_ice, 'V', -1. ) 528 529 #if defined key_agrif 599 530 CALL agrif_interp_lim3( 'U' ) 600 531 CALL agrif_interp_lim3( 'V' ) … … 605 536 #endif 606 537 607 DO jj = k_j1+1, k_jpj-1 538 !------------------------------------------------------------------------------! 539 ! 5) Recompute delta, shear and div (inputs for mechanical redistribution) 540 !------------------------------------------------------------------------------! 541 542 ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 543 DO jj = 2, jpjm1 608 544 DO ji = fs_2, fs_jpim1 609 IF ( vt_i(ji,jj) <= zvmin ) THEN 610 v_ice1(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji, jj-1) ) * e1t(ji+1,jj) & 611 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) & 612 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 613 614 u_ice2(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) & 615 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) & 616 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 617 ENDIF 545 546 ! divergence at T points 547 divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 548 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 549 & ) * r1_e12t(ji,jj) 550 551 ! tension at T points 552 zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & 553 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 554 & ) * r1_e12t(ji,jj) 555 556 ! shear at F points 557 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) & 558 & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 559 & ) * r1_e12f(ji,jj) 560 618 561 END DO 619 562 END DO 620 621 CALL lbc_lnk_multi( u_ice2(:,:), 'V', -1., v_ice1(:,:), 'U', -1. ) 622 623 ! Recompute delta, shear and div, inputs for mechanical redistribution 624 DO jj = k_j1+1, k_jpj-1 625 DO ji = fs_2, jpim1 !RB bug no vect opt due to zmask 626 !- divu_i(:,:), zdt(:,:): divergence and tension at centre 627 !- zds(:,:): shear on northeast corner of grid cells 628 IF ( vt_i(ji,jj) <= zvmin ) THEN 629 630 divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj ) * u_ice(ji-1,jj ) & 631 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji ,jj-1) * v_ice(ji ,jj-1) & 632 & ) * r1_e12t(ji,jj) 633 634 zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & 635 & -( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 636 & ) * r1_e12t(ji,jj) 637 ! 638 ! SB modif because ocean has no slip boundary condition 639 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) & 640 & +( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 641 & ) * r1_e12f(ji,jj) * ( 2.0 - fmask(ji,jj,1) ) & 642 & * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1) 643 644 zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj ) * v_ice1(ji-1,jj ) & 645 & + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji ,jj-1) * u_ice2(ji ,jj-1) ) * r1_e12t(ji,jj) 646 647 delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 ) 648 delta_i(ji,jj) = delta + rn_creepl 649 650 ENDIF 563 CALL lbc_lnk_multi( divu_i, 'T', 1., zds, 'F', 1. ) 564 565 566 DO jj = 2, jpjm1 567 DO ji = 2, jpim1 ! no vector loop 568 569 ! shear**2 at T points (doc eq. A16) 570 zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e12f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e12f(ji-1,jj ) & 571 & + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1) & 572 & ) * 0.25_wp * r1_e12t(ji,jj) 573 574 ! delta at T points 575 delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zds2 ) * usecc2 ) 576 delta_i(ji,jj) = delta + rn_creepl 577 578 shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zds2 ) 651 579 END DO 652 580 END DO 653 ! 654 !------------------------------------------------------------------------------! 655 ! 5) Store stress tensor and its invariants 656 !------------------------------------------------------------------------------! 657 ! * Invariants of the stress tensor are required for limitd_me 658 ! (accelerates convergence and improves stability) 659 DO jj = k_j1+1, k_jpj-1 660 DO ji = fs_2, fs_jpim1 661 zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u( ji-1, jj ) * v_ice1(ji-1,jj) & 662 & + e1v(ji,jj) * u_ice2(ji,jj) - e1v( ji , jj-1 ) * u_ice2(ji,jj-1) ) * r1_e12t(ji,jj) 663 shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst ) 664 END DO 665 END DO 666 667 ! Lateral boundary condition 668 CALL lbc_lnk_multi( divu_i (:,:), 'T', 1., delta_i(:,:), 'T', 1., shear_i(:,:), 'T', 1. ) 669 670 ! * Store the stress tensor for the next time step 581 CALL lbc_lnk_multi( delta_i, 'T', 1., shear_i, 'T', 1. ) 582 583 ! --- Store the stress tensor for the next time step --- ! 671 584 stress1_i (:,:) = zs1 (:,:) 672 585 stress2_i (:,:) = zs2 (:,:) 673 586 stress12_i(:,:) = zs12(:,:) 674 675 ! 587 ! 588 676 589 !------------------------------------------------------------------------------! 677 590 ! 6) Control prints of residual and charge ellipse … … 695 608 WRITE(charout,FMT="('lim_rhg :', I4, I6, I1, I1, A10)") 1000, numit, 0, 0, ' ch. ell. ' 696 609 CALL prt_ctl_info(charout) 697 DO jj = k_j1+1, k_jpj-1610 DO jj = 2, jpjm1 698 611 DO ji = 2, jpim1 699 612 IF (zpresh(ji,jj) > 1.0) THEN 700 sigma1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5) / ( 2*zpresh(ji,jj) )701 sigma2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5) / ( 2*zpresh(ji,jj) )613 zsig1 = ( zs1(ji,jj) + SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*zpresh(ji,jj) ) 614 zsig2 = ( zs1(ji,jj) - SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*zpresh(ji,jj) ) 702 615 WRITE(charout,FMT="('lim_rhg :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)") 703 616 CALL prt_ctl_info(charout) … … 710 623 ENDIF 711 624 ! 712 CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 713 CALL wrk_dealloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask ) 714 CALL wrk_dealloc( jpi,jpj, zf1 , zu_ice, zf2 , zv_ice , zdt , zds ) 715 CALL wrk_dealloc( jpi,jpj, zs1 , zs2 , zs12 , zresr , zpice ) 625 CALL wrk_dealloc( jpi,jpj, zpresh, z1_e1t0, z1_e2t0, zp_delt ) 626 CALL wrk_dealloc( jpi,jpj, za1, za2, zmass1, zmass2, zcor1, zcor2 ) 627 CALL wrk_dealloc( jpi,jpj, zspgu, zspgv, v_oce1, u_oce2, v_ice1, u_ice2, zf1, zf2 ) 628 CALL wrk_dealloc( jpi,jpj, zds, zdt, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 629 CALL wrk_dealloc( jpi,jpj, zswitch1, zswitch2 ) 716 630 717 631 END SUBROUTINE lim_rhg … … 722 636 !!---------------------------------------------------------------------- 723 637 CONTAINS 724 SUBROUTINE lim_rhg ( k1 , k2 )! Dummy routine725 WRITE(*,*) 'lim_rhg: You should not have seen this print! error?' , k1, k2638 SUBROUTINE lim_rhg ! Dummy routine 639 WRITE(*,*) 'lim_rhg: You should not have seen this print! error?' 726 640 END SUBROUTINE lim_rhg 727 641 #endif
Note: See TracChangeset
for help on using the changeset viewer.
