- Timestamp:
- 2016-01-08T10:35:19+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
r4688 r6225 50 50 PUBLIC lim_rhg ! routine called by lim_dyn (or lim_dyn_2) 51 51 52 REAL(wp) :: epsi10 = 1.e-10_wp !53 54 52 !! * Substitutions 55 53 # include "vectopt_loop_substitute.h90" … … 104 102 !! and charge ellipse. 105 103 !! The user should make sure that the parameters 106 !! n evp, telast andcreepl maintain stress state104 !! nn_nevp, elastic time scale and rn_creepl maintain stress state 107 105 !! on the charge ellipse for plastic flow 108 106 !! e.g. in the Canadian Archipelago … … 110 108 !! References : Hunke and Dukowicz, JPO97 111 109 !! Bouillon et al., Ocean Modelling 2009 112 !! Vancoppenolle et al., Ocean Modelling 2008113 110 !!------------------------------------------------------------------- 114 111 INTEGER, INTENT(in) :: k_j1 ! southern j-index for ice computation … … 119 116 CHARACTER (len=50) :: charout 120 117 REAL(wp) :: zt11, zt12, zt21, zt22, ztagnx, ztagny, delta ! 121 REAL(wp) :: za, zstms, zsang, zmask ! local scalars 122 123 REAL(wp) :: dtevp ! time step for subcycling 124 REAL(wp) :: dtotel, ecc2, ecci ! square of yield ellipse eccenticity 125 REAL(wp) :: z0, zr, zcca, zccb ! temporary scalars 126 REAL(wp) :: zu_ice2, zv_ice1 ! 127 REAL(wp) :: zddc, zdtc, zzdst ! delta on corners and on centre 128 REAL(wp) :: zdsshx, zdsshy ! term for the gradient of ocean surface 129 REAL(wp) :: sigma1, sigma2 ! internal ice stress 118 REAL(wp) :: za, zstms ! local scalars 119 REAL(wp) :: zc1, zc2, zc3 ! ice mass 120 121 REAL(wp) :: dtevp , z1_dtevp ! time step for subcycling 122 REAL(wp) :: dtotel, z1_dtotel, ecc2, ecci ! square of yield ellipse eccenticity 123 REAL(wp) :: z0, zr, zcca, zccb ! temporary scalars 124 REAL(wp) :: zu_ice2, zv_ice1 ! 125 REAL(wp) :: zddc, zdtc ! delta on corners and on centre 126 REAL(wp) :: zdst ! shear at the center of the grid point 127 REAL(wp) :: zdsshx, zdsshy ! term for the gradient of ocean surface 128 REAL(wp) :: sigma1, sigma2 ! internal ice stress 130 129 131 130 REAL(wp) :: zresm ! Maximal error on ice velocity 132 REAL(wp) :: zindb ! ice (1) or not (0)133 REAL(wp) :: zdummy ! dummy argument134 131 REAL(wp) :: zintb, zintn ! dummy argument 135 132 … … 140 137 REAL(wp), POINTER, DIMENSION(:,:) :: zcorl1, zcorl2 ! coriolis parameter on U/V points 141 138 REAL(wp), POINTER, DIMENSION(:,:) :: za1ct , za2ct ! temporary arrays 142 REAL(wp), POINTER, DIMENSION(:,:) :: zc1 ! ice mass 143 REAL(wp), POINTER, DIMENSION(:,:) :: zusw ! temporary weight for ice strength computation 144 REAL(wp), POINTER, DIMENSION(:,:) :: u_oce1, v_oce1 ! ocean u/v component on U points 145 REAL(wp), POINTER, DIMENSION(:,:) :: u_oce2, v_oce2 ! ocean u/v component on V points 139 REAL(wp), POINTER, DIMENSION(:,:) :: v_oce1 ! ocean u/v component on U points 140 REAL(wp), POINTER, DIMENSION(:,:) :: u_oce2 ! ocean u/v component on V points 146 141 REAL(wp), POINTER, DIMENSION(:,:) :: u_ice2, v_ice1 ! ice u/v component on V/U point 147 142 REAL(wp), POINTER, DIMENSION(:,:) :: zf1 , zf2 ! arrays for internal stresses 143 REAL(wp), POINTER, DIMENSION(:,:) :: zmask ! mask ocean grid points 148 144 149 REAL(wp), POINTER, DIMENSION(:,:) :: zd d , zdt ! Divergence andtension at centre of grid cells145 REAL(wp), POINTER, DIMENSION(:,:) :: zdt ! tension at centre of grid cells 150 146 REAL(wp), POINTER, DIMENSION(:,:) :: zds ! Shear on northeast corner of grid cells 151 REAL(wp), POINTER, DIMENSION(:,:) :: zdst ! Shear on centre of grid cells152 REAL(wp), POINTER, DIMENSION(:,:) :: deltat, deltac ! Delta at centre and corners of grid cells153 147 REAL(wp), POINTER, DIMENSION(:,:) :: zs1 , zs2 ! Diagonal stress tensor components zs1 and zs2 154 148 REAL(wp), POINTER, DIMENSION(:,:) :: zs12 ! Non-diagonal stress tensor component zs12 … … 157 151 ! ocean surface (ssh_m) if ice is not embedded 158 152 ! ice top surface if ice is embedded 153 154 REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter 155 REAL(wp), PARAMETER :: zvmin = 1.0e-03_wp ! ice volume below which ice velocity equals ocean velocity 159 156 !!------------------------------------------------------------------- 160 157 161 158 CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 162 CALL wrk_alloc( jpi,jpj, zc1 , u_oce1, u_oce2, u_ice2, zusw , v_oce1 , v_oce2, v_ice1)163 CALL wrk_alloc( jpi,jpj, zf1 , deltat, zu_ice, zf2 , deltac, zv_ice , zdd , zdt , zds , zdst)164 CALL wrk_alloc( jpi,jpj, zd d , zdt , zds , zs1 , zs2 , zs12 , zresr , zpice )159 CALL wrk_alloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask ) 160 CALL wrk_alloc( jpi,jpj, zf1 , zu_ice, zf2 , zv_ice , zdt , zds ) 161 CALL wrk_alloc( jpi,jpj, zdt , zds , zs1 , zs2 , zs12 , zresr , zpice ) 165 162 166 163 #if defined key_lim2 && ! defined key_lim2_vp 167 164 # if defined key_agrif 168 USE ice_2, vt_s => hsnm169 USE ice_2, vt_i => hicm165 USE ice_2, vt_s => hsnm 166 USE ice_2, vt_i => hicm 170 167 # else 171 vt_s => hsnm172 vt_i => hicm168 vt_s => hsnm 169 vt_i => hicm 173 170 # endif 174 at_i(:,:) = 1. - frld(:,:)171 at_i(:,:) = 1. - frld(:,:) 175 172 #endif 176 173 #if defined key_agrif && defined key_lim2 177 CALL agrif_rhg_lim2_load ! First interpolation of coarse values178 #endif 179 ! 180 !------------------------------------------------------------------------------! 181 ! 1) Ice -Snow mass (zc1), icestrength (zpresh) !174 CALL agrif_rhg_lim2_load ! First interpolation of coarse values 175 #endif 176 ! 177 !------------------------------------------------------------------------------! 178 ! 1) Ice strength (zpresh) ! 182 179 !------------------------------------------------------------------------------! 183 180 ! 184 181 ! Put every vector to 0 185 zpresh (:,:) = 0._wp ; zc1 (:,:) = 0._wp 182 delta_i(:,:) = 0._wp ; 183 zpresh (:,:) = 0._wp ; 186 184 zpreshc(:,:) = 0._wp 187 185 u_ice2 (:,:) = 0._wp ; v_ice1(:,:) = 0._wp 188 zdd (:,:) = 0._wp ; zdt (:,:) = 0._wp ; zds(:,:) = 0._wp 186 divu_i (:,:) = 0._wp ; zdt (:,:) = 0._wp ; zds(:,:) = 0._wp 187 shear_i(:,:) = 0._wp 189 188 190 189 #if defined key_lim3 191 CALL lim_itd_me_icestrength( ridge_scheme_swi ) ! LIM-3: Ice strength on T-points 192 #endif 193 194 !CDIR NOVERRCHK 190 CALL lim_itd_me_icestrength( nn_icestr ) ! LIM-3: Ice strength on T-points 191 #endif 192 195 193 DO jj = k_j1 , k_jpj ! Ice mass and temp variables 196 !CDIR NOVERRCHK197 194 DO ji = 1 , jpi 198 zc1(ji,jj) = tms(ji,jj) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj) )199 195 #if defined key_lim3 200 zpresh(ji,jj) = tm s(ji,jj) * strength(ji,jj)196 zpresh(ji,jj) = tmask(ji,jj,1) * strength(ji,jj) 201 197 #endif 202 198 #if defined key_lim2 203 zpresh(ji,jj) = tm s(ji,jj) * pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) )204 #endif 205 ! tmi= 1 where there is ice or on land206 tmi(ji,jj) = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - epsd ) ) ) * tms(ji,jj)199 zpresh(ji,jj) = tmask(ji,jj,1) * pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) ) 200 #endif 201 ! zmask = 1 where there is ice or on land 202 zmask(ji,jj) = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - zepsi ) ) ) * tmask(ji,jj,1) 207 203 END DO 208 204 END DO … … 210 206 ! Ice strength on grid cell corners (zpreshc) 211 207 ! needed for calculation of shear stress 212 !CDIR NOVERRCHK213 208 DO jj = k_j1+1, k_jpj-1 214 !CDIR NOVERRCHK215 209 DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1) 216 zstms = tms(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + & 217 & tms(ji,jj+1) * wght(ji+1,jj+1,1,2) + & 218 & tms(ji+1,jj) * wght(ji+1,jj+1,2,1) + & 219 & tms(ji,jj) * wght(ji+1,jj+1,1,1) 220 zusw(ji,jj) = 1.0 / MAX( zstms, epsd ) 221 zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + & 222 & zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) + & 223 & zpresh(ji+1,jj) * wght(ji+1,jj+1,2,1) + & 224 & zpresh(ji,jj) * wght(ji+1,jj+1,1,1) & 225 & ) * zusw(ji,jj) 210 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) + & 211 & tmask(ji+1,jj,1) * wght(ji+1,jj+1,2,1) + tmask(ji,jj,1) * wght(ji+1,jj+1,1,1) 212 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) + & 213 & zpresh(ji+1,jj) * wght(ji+1,jj+1,2,1) + zpresh(ji,jj) * wght(ji+1,jj+1,1,1) & 214 & ) / MAX( zstms, zepsi ) 226 215 END DO 227 216 END DO 228 229 217 CALL lbc_lnk( zpreshc(:,:), 'F', 1. ) 230 218 ! … … 241 229 ! zcorl2: Coriolis parameter on V-points 242 230 ! (ztagnx,ztagny): wind stress on U/V points 243 ! u_oce1: ocean u component on u points244 231 ! v_oce1: ocean v component on u points 245 232 ! u_oce2: ocean u component on v points 246 ! v_oce2: ocean v component on v points247 233 248 234 IF( nn_ice_embd == 2 ) THEN !== embedded sea ice: compute representative ice top surface ==! 249 250 251 235 ! 236 ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1} 237 ! = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1} 252 238 zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp 253 254 255 239 ! 240 ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1} 241 ! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1}) 256 242 zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 257 243 ! 258 244 zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0 259 245 ! 260 246 ELSE !== non-embedded sea ice: use ocean surface for slope calculation ==! 261 247 zpice(:,:) = ssh_m(:,:) … … 265 251 DO ji = fs_2, fs_jpim1 266 252 267 zt11 = tms(ji ,jj) * e1t(ji ,jj) 268 zt12 = tms(ji+1,jj) * e1t(ji+1,jj) 269 zt21 = tms(ji,jj ) * e2t(ji,jj ) 270 zt22 = tms(ji,jj+1) * e2t(ji,jj+1) 253 zc1 = tmask(ji ,jj ,1) * ( rhosn * vt_s(ji ,jj ) + rhoic * vt_i(ji ,jj ) ) 254 zc2 = tmask(ji+1,jj ,1) * ( rhosn * vt_s(ji+1,jj ) + rhoic * vt_i(ji+1,jj ) ) 255 zc3 = tmask(ji ,jj+1,1) * ( rhosn * vt_s(ji ,jj+1) + rhoic * vt_i(ji ,jj+1) ) 256 257 zt11 = tmask(ji ,jj,1) * e1t(ji ,jj) 258 zt12 = tmask(ji+1,jj,1) * e1t(ji+1,jj) 259 zt21 = tmask(ji,jj ,1) * e2t(ji,jj ) 260 zt22 = tmask(ji,jj+1,1) * e2t(ji,jj+1) 271 261 272 262 ! Leads area. 273 zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + epsd)274 zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + epsd)263 zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + zepsi ) 264 zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + zepsi ) 275 265 276 266 ! Mass, coriolis coeff. and currents 277 zmass1(ji,jj) = ( zt12 *zc1(ji,jj) + zt11*zc1(ji+1,jj) ) / (zt11+zt12+epsd)278 zmass2(ji,jj) = ( zt22 *zc1(ji,jj) + zt21*zc1(ji,jj+1) ) / (zt21+zt22+epsd)279 zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj) *fcor(ji,jj) + e1t(ji,jj)*fcor(ji+1,jj) ) &280 & / ( e1t(ji,jj) + e1t(ji+1,jj) + epsd)281 zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1) *fcor(ji,jj) + e2t(ji,jj)*fcor(ji,jj+1) ) &282 & / ( e2t(ji,jj+1) + e2t(ji,jj) + epsd)267 zmass1(ji,jj) = ( zt12 * zc1 + zt11 * zc2 ) / ( zt11 + zt12 + zepsi ) 268 zmass2(ji,jj) = ( zt22 * zc1 + zt21 * zc3 ) / ( zt21 + zt22 + zepsi ) 269 zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj) * fcor(ji,jj) + e1t(ji,jj) * fcor(ji+1,jj) ) & 270 & / ( e1t(ji,jj) + e1t(ji+1,jj) + zepsi ) 271 zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1) * fcor(ji,jj) + e2t(ji,jj) * fcor(ji,jj+1) ) & 272 & / ( e2t(ji,jj+1) + e2t(ji,jj) + zepsi ) 283 273 ! 284 u_oce1(ji,jj) = u_oce(ji,jj)285 v_oce2(ji,jj) = v_oce(ji,jj)286 287 274 ! Ocean has no slip boundary condition 288 v_oce1(ji,jj) = 0.5 *( (v_oce(ji,jj)+v_oce(ji,jj-1))*e1t(ji,jj)&289 & +(v_oce(ji+1,jj)+v_oce(ji+1,jj-1))*e1t(ji+1,jj))&290 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)291 292 u_oce2(ji,jj) = 0.5 *((u_oce(ji,jj)+u_oce(ji-1,jj))*e2t(ji,jj)&293 & +(u_oce(ji,jj+1)+u_oce(ji-1,jj+1))*e2t(ji,jj+1))&294 & / (e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)275 v_oce1(ji,jj) = 0.5 * ( ( v_oce(ji ,jj) + v_oce(ji ,jj-1) ) * e1t(ji,jj) & 276 & + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji+1,jj) ) & 277 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 278 279 u_oce2(ji,jj) = 0.5 * ( ( u_oce(ji,jj ) + u_oce(ji-1,jj ) ) * e2t(ji,jj) & 280 & + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj+1) ) & 281 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 295 282 296 283 ! Wind stress at U,V-point … … 304 291 ! include it later 305 292 306 zdsshx = ( zpice(ji+1,jj) - zpice(ji,jj) ) /e1u(ji,jj)307 zdsshy = ( zpice(ji,jj+1) - zpice(ji,jj) ) /e2v(ji,jj)293 zdsshx = ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj) 294 zdsshy = ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj) 308 295 309 296 za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx … … 319 306 ! 320 307 ! Time step for subcycling 321 dtevp = rdt_ice / nevp 308 dtevp = rdt_ice / nn_nevp 309 #if defined key_lim3 310 dtotel = dtevp / ( 2._wp * rn_relast * rdt_ice ) 311 #else 322 312 dtotel = dtevp / ( 2._wp * telast ) 323 313 #endif 314 z1_dtotel = 1._wp / ( 1._wp + dtotel ) 315 z1_dtevp = 1._wp / dtevp 324 316 !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter) 325 ecc2 = ecc *ecc317 ecc2 = rn_ecc * rn_ecc 326 318 ecci = 1. / ecc2 327 319 … … 332 324 333 325 ! !----------------------! 334 DO jter = 1 , n evp! loop over jter !326 DO jter = 1 , nn_nevp ! loop over jter ! 335 327 ! !----------------------! 336 328 DO jj = k_j1, k_jpj-1 … … 340 332 341 333 DO jj = k_j1+1, k_jpj-1 342 DO ji = fs_2, jpim1 !RB bug no vect opt due to tmi334 DO ji = fs_2, fs_jpim1 !RB bug no vect opt due to zmask 343 335 344 336 ! 345 337 !- Divergence, tension and shear (Section a. Appendix B of Hunke & Dukowicz, 2002) 346 !- zdd(:,:), zdt(:,:): divergence and tension at centre of grid cells338 !- divu_i(:,:), zdt(:,:): divergence and tension at centre of grid cells 347 339 !- zds(:,:): shear on northeast corner of grid cells 348 340 ! … … 353 345 ! bugs (Martin, for Miguel). 354 346 ! 355 !- ALSO: arrays zd d, zdt, zds and delta could347 !- ALSO: arrays zdt, zds and delta could 356 348 ! be removed in the future to minimise memory demand. 357 349 ! … … 361 353 ! 362 354 ! 363 zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) & 364 & -e2u(ji-1,jj)*u_ice(ji-1,jj) & 365 & +e1v(ji,jj)*v_ice(ji,jj) & 366 & -e1v(ji,jj-1)*v_ice(ji,jj-1) & 367 & ) & 368 & / area(ji,jj) 369 370 zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj) & 371 & -u_ice(ji-1,jj)/e2u(ji-1,jj) & 372 & )*e2t(ji,jj)*e2t(ji,jj) & 373 & -( v_ice(ji,jj)/e1v(ji,jj) & 374 & -v_ice(ji,jj-1)/e1v(ji,jj-1) & 375 & )*e1t(ji,jj)*e1t(ji,jj) & 376 & ) & 377 & / area(ji,jj) 355 divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 356 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 357 & ) * r1_e1e2t(ji,jj) 358 359 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) & 360 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 361 & ) * r1_e1e2t(ji,jj) 378 362 379 363 ! 380 zds(ji,jj) = ( ( u_ice(ji,jj+1)/e1u(ji,jj+1) & 381 & -u_ice(ji,jj)/e1u(ji,jj) & 382 & )*e1f(ji,jj)*e1f(ji,jj) & 383 & +( v_ice(ji+1,jj)/e2v(ji+1,jj) & 384 & -v_ice(ji,jj)/e2v(ji,jj) & 385 & )*e2f(ji,jj)*e2f(ji,jj) & 386 & ) & 387 & / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) & 388 & * tmi(ji,jj) * tmi(ji,jj+1) & 389 & * tmi(ji+1,jj) * tmi(ji+1,jj+1) 390 391 392 v_ice1(ji,jj) = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj) & 393 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & 394 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 395 396 u_ice2(ji,jj) = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1) & 397 & +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) & 398 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 399 400 END DO 401 END DO 402 CALL lbc_lnk( v_ice1, 'U', -1. ) ; CALL lbc_lnk( u_ice2, 'V', -1. ) ! lateral boundary cond. 403 404 !CDIR NOVERRCHK 364 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) & 365 & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 366 & ) * r1_e1e2f(ji,jj) * ( 2._wp - fmask(ji,jj,1) ) & 367 & * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1) 368 369 370 v_ice1(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji+1,jj) & 371 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) & 372 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 373 374 u_ice2(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) & 375 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) & 376 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 377 END DO 378 END DO 379 380 CALL lbc_lnk_multi( v_ice1, 'U', -1., u_ice2, 'V', -1. ) ! lateral boundary cond. 381 405 382 DO jj = k_j1+1, k_jpj-1 406 !CDIR NOVERRCHK407 383 DO ji = fs_2, fs_jpim1 408 384 409 385 !- Calculate Delta at centre of grid cells 410 zzdst = ( e2u(ji , jj) * v_ice1(ji ,jj) & 411 & - e2u(ji-1, jj) * v_ice1(ji-1,jj) & 412 & + e1v(ji, jj ) * u_ice2(ji,jj ) & 413 & - e1v(ji, jj-1) * u_ice2(ji,jj-1) & 414 & ) & 415 & / area(ji,jj) 416 417 delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zzdst*zzdst ) * usecc2 ) 418 ! MV rewriting 419 ! deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 + (zdt(ji,jj)**2 + zzdst**2)*usecc2), creepl ) 420 !!gm faster to replace the line above with simply: 421 !! deltat(ji,jj) = MAX( delta, creepl ) 422 !!gm end 423 deltat(ji,jj) = delta + creepl 424 ! END MV 425 !-Calculate stress tensor components zs1 and zs2 426 !-at centre of grid cells (see section 3.5 of CICE user's guide). 427 !zs1(ji,jj) = ( zs1(ji,jj) - dtotel*( ( 1._wp - alphaevp) * zs1(ji,jj) + & 428 ! & ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) * zpresh(ji,jj) ) ) & 429 ! & / ( 1._wp + alphaevp * dtotel ) 430 431 !zs2(ji,jj) = ( zs2(ji,jj) - dtotel * ( ( 1._wp - alphaevp ) * ecc2 * zs2(ji,jj) - & 432 ! zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) ) & 433 ! & / ( 1._wp + alphaevp * ecc2 * dtotel ) 434 435 ! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp) 436 zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( zdd(ji,jj) / deltat(ji,jj) - delta / deltat(ji,jj) ) & 437 & * zpresh(ji,jj) ) ) / ( 1._wp + dtotel ) 438 zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) ) & 439 & / ( 1._wp + dtotel ) 440 441 END DO 442 END DO 443 444 CALL lbc_lnk( zs1(:,:), 'T', 1. ) 445 CALL lbc_lnk( zs2(:,:), 'T', 1. ) 446 447 !CDIR NOVERRCHK 448 DO jj = k_j1+1, k_jpj-1 449 !CDIR NOVERRCHK 450 DO ji = fs_2, fs_jpim1 386 zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj ) * v_ice1(ji-1,jj ) & 387 & + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji ,jj-1) * u_ice2(ji ,jj-1) & 388 & ) * r1_e1e2t(ji,jj) 389 390 delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 ) 391 delta_i(ji,jj) = delta + rn_creepl 392 451 393 !- Calculate Delta on corners 452 zddc = ( ( v_ice1(ji,jj+1)/e1u(ji,jj+1) & 453 & -v_ice1(ji,jj)/e1u(ji,jj) & 454 & )*e1f(ji,jj)*e1f(ji,jj) & 455 & +( u_ice2(ji+1,jj)/e2v(ji+1,jj) & 456 & -u_ice2(ji,jj)/e2v(ji,jj) & 457 & )*e2f(ji,jj)*e2f(ji,jj) & 458 & ) & 459 & / ( e1f(ji,jj) * e2f(ji,jj) ) 460 461 zdtc = (-( v_ice1(ji,jj+1)/e1u(ji,jj+1) & 462 & -v_ice1(ji,jj)/e1u(ji,jj) & 463 & )*e1f(ji,jj)*e1f(ji,jj) & 464 & +( u_ice2(ji+1,jj)/e2v(ji+1,jj) & 465 & -u_ice2(ji,jj)/e2v(ji,jj) & 466 & )*e2f(ji,jj)*e2f(ji,jj) & 467 & ) & 468 & / ( e1f(ji,jj) * e2f(ji,jj) ) 469 470 deltac(ji,jj) = SQRT(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl 471 472 !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide). 473 !zs12(ji,jj) = ( zs12(ji,jj) - dtotel * ( (1.0-alphaevp) * ecc2 * zs12(ji,jj) - zds(ji,jj) / & 474 ! & ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) ) & 475 ! & / ( 1._wp + alphaevp * ecc2 * dtotel ) 476 477 ! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp) 478 zs12(ji,jj) = ( zs12(ji,jj) + dtotel * & 479 & ( ecci * zds(ji,jj) / ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) ) & 480 & / ( 1.0 + dtotel ) 481 482 END DO ! ji 483 END DO ! jj 484 485 CALL lbc_lnk( zs12(:,:), 'F', 1. ) 486 394 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) & 395 & + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 396 & ) * r1_e1e2f(ji,jj) 397 398 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) & 399 & + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 400 & ) * r1_e1e2f(ji,jj) 401 402 zddc = SQRT( zddc**2 + ( zdtc**2 + zds(ji,jj)**2 ) * usecc2 ) + rn_creepl 403 404 !-Calculate stress tensor components zs1 and zs2 at centre of grid cells (see section 3.5 of CICE user's guide). 405 zs1(ji,jj) = ( zs1 (ji,jj) + dtotel * ( divu_i(ji,jj) - delta ) / delta_i(ji,jj) * zpresh(ji,jj) & 406 & ) * z1_dtotel 407 zs2(ji,jj) = ( zs2 (ji,jj) + dtotel * ecci * zdt(ji,jj) / delta_i(ji,jj) * zpresh(ji,jj) & 408 & ) * z1_dtotel 409 !-Calculate stress tensor component zs12 at corners 410 zs12(ji,jj) = ( zs12(ji,jj) + dtotel * ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj) & 411 & ) * z1_dtotel 412 413 END DO 414 END DO 415 416 CALL lbc_lnk_multi( zs1 , 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 417 487 418 ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) 488 419 DO jj = k_j1+1, k_jpj-1 489 420 DO ji = fs_2, fs_jpim1 490 421 !- contribution of zs1, zs2 and zs12 to zf1 491 zf1(ji,jj) = 0.5 *( (zs1(ji+1,jj)-zs1(ji,jj))*e2u(ji,jj)&492 & +(zs2(ji+1,jj)*e2t(ji+1,jj)**2-zs2(ji,jj)*e2t(ji,jj)**2)/e2u(ji,jj)&493 & +2.0*(zs12(ji,jj)*e1f(ji,jj)**2-zs12(ji,jj-1)*e1f(ji,jj-1)**2)/e1u(ji,jj)&494 & ) / ( e1u(ji,jj)*e2u(ji,jj))422 zf1(ji,jj) = 0.5 * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) & 423 & + ( zs2(ji+1,jj) * e2t(ji+1,jj)**2 - zs2(ji,jj) * e2t(ji,jj)**2 ) * r1_e2u(ji,jj) & 424 & + 2.0 * ( zs12(ji,jj) * e1f(ji,jj)**2 - zs12(ji,jj-1) * e1f(ji,jj-1)**2 ) * r1_e1u(ji,jj) & 425 & ) * r1_e1e2u(ji,jj) 495 426 ! contribution of zs1, zs2 and zs12 to zf2 496 zf2(ji,jj) = 0.5*( (zs1(ji,jj+1)-zs1(ji,jj))*e1v(ji,jj) & 497 & -(zs2(ji,jj+1)*e1t(ji,jj+1)**2 - zs2(ji,jj)*e1t(ji,jj)**2)/e1v(ji,jj) & 498 & + 2.0*(zs12(ji,jj)*e2f(ji,jj)**2 - & 499 zs12(ji-1,jj)*e2f(ji-1,jj)**2)/e2v(ji,jj) & 500 & ) / ( e1v(ji,jj)*e2v(ji,jj) ) 427 zf2(ji,jj) = 0.5 * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) & 428 & - ( zs2(ji,jj+1) * e1t(ji,jj+1)**2 - zs2(ji,jj) * e1t(ji,jj)**2 ) * r1_e1v(ji,jj) & 429 & + 2.0 * ( zs12(ji,jj) * e2f(ji,jj)**2 - zs12(ji-1,jj) * e2f(ji-1,jj)**2 ) * r1_e2v(ji,jj) & 430 & ) * r1_e1e2v(ji,jj) 501 431 END DO 502 432 END DO … … 508 438 IF (MOD(jter,2).eq.0) THEN 509 439 510 !CDIR NOVERRCHK511 440 DO jj = k_j1+1, k_jpj-1 512 !CDIR NOVERRCHK513 441 DO ji = fs_2, fs_jpim1 514 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj) 515 zsang = SIGN ( 1.0 , fcor(ji,jj) ) * sangvg 516 z0 = zmass1(ji,jj)/dtevp 442 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1) 443 z0 = zmass1(ji,jj) * z1_dtevp 517 444 518 445 ! SB modif because ocean has no slip boundary condition 519 zv_ice1 = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj) & 520 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj)) & 521 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 522 za = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + & 523 (zv_ice1-v_oce1(ji,jj))**2) * (1.0-zfrld1(ji,jj)) 524 zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + & 525 za*(cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj)) 526 zcca = z0+za*cangvg 527 zccb = zcorl1(ji,jj)+za*zsang 528 u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask 529 446 zv_ice1 = 0.5 * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji ,jj) & 447 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) ) & 448 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 449 za = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 + & 450 & ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) ) 451 zr = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj) 452 zcca = z0 + za 453 zccb = zcorl1(ji,jj) 454 u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch 530 455 END DO 531 456 END DO … … 533 458 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 534 459 #if defined key_agrif && defined key_lim2 535 CALL agrif_rhg_lim2( jter, n evp, 'U' )460 CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 536 461 #endif 537 462 #if defined key_bdy 538 ! clem: change u_ice and v_ice at the boundary for each iteration539 463 CALL bdy_ice_lim_dyn( 'U' ) 540 464 #endif 541 465 542 !CDIR NOVERRCHK543 466 DO jj = k_j1+1, k_jpj-1 544 !CDIR NOVERRCHK545 467 DO ji = fs_2, fs_jpim1 546 468 547 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj) 548 zsang = SIGN(1.0,fcor(ji,jj))*sangvg 549 z0 = zmass2(ji,jj)/dtevp 469 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1) 470 z0 = zmass2(ji,jj) * z1_dtevp 550 471 ! SB modif because ocean has no slip boundary condition 551 zu_ice2 = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj) & 552 & + (u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1)) & 553 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 554 za = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + & 555 (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 556 zr = z0*v_ice(ji,jj) + zf2(ji,jj) + & 557 za2ct(ji,jj) + za*(cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj)) 558 zcca = z0+za*cangvg 559 zccb = zcorl2(ji,jj)+za*zsang 560 v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask 561 472 zu_ice2 = 0.5 * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj) & 473 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) ) & 474 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 475 za = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 + & 476 & ( v_ice(ji,jj) - v_oce(ji,jj))**2 ) * ( 1.0 - zfrld2(ji,jj) ) 477 zr = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj) 478 zcca = z0 + za 479 zccb = zcorl2(ji,jj) 480 v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch 562 481 END DO 563 482 END DO … … 565 484 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 566 485 #if defined key_agrif && defined key_lim2 567 CALL agrif_rhg_lim2( jter, n evp, 'V' )486 CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 568 487 #endif 569 488 #if defined key_bdy 570 ! clem: change u_ice and v_ice at the boundary for each iteration571 489 CALL bdy_ice_lim_dyn( 'V' ) 572 490 #endif 573 491 574 492 ELSE 575 !CDIR NOVERRCHK576 493 DO jj = k_j1+1, k_jpj-1 577 !CDIR NOVERRCHK578 494 DO ji = fs_2, fs_jpim1 579 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj) 580 zsang = SIGN(1.0,fcor(ji,jj))*sangvg 581 z0 = zmass2(ji,jj)/dtevp 495 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1) 496 z0 = zmass2(ji,jj) * z1_dtevp 582 497 ! SB modif because ocean has no slip boundary condition 583 zu_ice2 = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj) & 584 & +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1)) & 585 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 586 587 za = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + & 588 (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 589 zr = z0*v_ice(ji,jj) + zf2(ji,jj) + & 590 za2ct(ji,jj) + za*(cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj)) 591 zcca = z0+za*cangvg 592 zccb = zcorl2(ji,jj)+za*zsang 593 v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask 594 498 zu_ice2 = 0.5 * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj) & 499 & +( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) ) & 500 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 501 502 za = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 + & 503 & ( v_ice(ji,jj) - v_oce(ji,jj) )**2 ) * ( 1.0 - zfrld2(ji,jj) ) 504 zr = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj) 505 zcca = z0 + za 506 zccb = zcorl2(ji,jj) 507 v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch 595 508 END DO 596 509 END DO … … 598 511 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 599 512 #if defined key_agrif && defined key_lim2 600 CALL agrif_rhg_lim2( jter, n evp, 'V' )513 CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 601 514 #endif 602 515 #if defined key_bdy 603 ! clem: change u_ice and v_ice at the boundary for each iteration604 516 CALL bdy_ice_lim_dyn( 'V' ) 605 517 #endif 606 518 607 !CDIR NOVERRCHK608 519 DO jj = k_j1+1, k_jpj-1 609 !CDIR NOVERRCHK610 520 DO ji = fs_2, fs_jpim1 611 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj) 612 zsang = SIGN(1.0,fcor(ji,jj))*sangvg 613 z0 = zmass1(ji,jj)/dtevp 614 ! SB modif because ocean has no slip boundary condition 615 ! GG Bug 616 ! zv_ice1 = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj) & 617 ! & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & 618 ! & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 619 zv_ice1 = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj) & 620 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj)) & 621 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 622 623 za = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + & 624 (zv_ice1-v_oce1(ji,jj))**2)*(1.0-zfrld1(ji,jj)) 625 zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + & 626 za*(cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj)) 627 zcca = z0+za*cangvg 628 zccb = zcorl1(ji,jj)+za*zsang 629 u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask 630 END DO ! ji 631 END DO ! jj 521 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1) 522 z0 = zmass1(ji,jj) * z1_dtevp 523 zv_ice1 = 0.5 * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji,jj) & 524 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) ) & 525 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 526 527 za = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 + & 528 & ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) ) 529 zr = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj) 530 zcca = z0 + za 531 zccb = zcorl1(ji,jj) 532 u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch 533 END DO 534 END DO 632 535 633 536 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 634 537 #if defined key_agrif && defined key_lim2 635 CALL agrif_rhg_lim2( jter, n evp, 'U' )538 CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 636 539 #endif 637 540 #if defined key_bdy 638 ! clem: change u_ice and v_ice at the boundary for each iteration639 541 CALL bdy_ice_lim_dyn( 'U' ) 640 542 #endif … … 645 547 !--- Convergence test. 646 548 DO jj = k_j1+1 , k_jpj-1 647 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) , & 648 ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 649 END DO 650 zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) ) 549 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 550 END DO 551 zresm = MAXVAL( zresr( 1:jpi, k_j1+1:k_jpj-1 ) ) 651 552 IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain 652 553 ENDIF … … 659 560 ! 4) Prevent ice velocities when the ice is thin 660 561 !------------------------------------------------------------------------------! 661 ! If the ice thickness is below hminrhg (5cm) then ice velocity should equal the 662 ! ocean velocity, 663 ! This prevents high velocity when ice is thin 664 !CDIR NOVERRCHK 562 ! If the ice volume is below zvmin then ice velocity should equal the 563 ! ocean velocity. This prevents high velocity when ice is thin 665 564 DO jj = k_j1+1, k_jpj-1 666 !CDIR NOVERRCHK667 565 DO ji = fs_2, fs_jpim1 668 zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) ) 669 !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 ) 670 zdummy = vt_i(ji,jj) 671 IF ( zdummy .LE. hminrhg ) THEN 566 IF ( vt_i(ji,jj) <= zvmin ) THEN 672 567 u_ice(ji,jj) = u_oce(ji,jj) 673 568 v_ice(ji,jj) = v_oce(ji,jj) 674 ENDIF ! zdummy569 ENDIF 675 570 END DO 676 571 END DO 677 572 678 CALL lbc_lnk ( u_ice(:,:), 'U', -1. )679 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 573 CALL lbc_lnk_multi( u_ice(:,:), 'U', -1., v_ice(:,:), 'V', -1. ) 574 680 575 #if defined key_agrif && defined key_lim2 681 CALL agrif_rhg_lim2( n evp ,nevp, 'U' )682 CALL agrif_rhg_lim2( n evp ,nevp, 'V' )576 CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'U' ) 577 CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'V' ) 683 578 #endif 684 579 #if defined key_bdy 685 ! clem: change u_ice and v_ice at the boundary686 580 CALL bdy_ice_lim_dyn( 'U' ) 687 581 CALL bdy_ice_lim_dyn( 'V' ) … … 690 584 DO jj = k_j1+1, k_jpj-1 691 585 DO ji = fs_2, fs_jpim1 692 zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) ) 693 !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 ) 694 zdummy = vt_i(ji,jj) 695 IF ( zdummy .LE. hminrhg ) THEN 696 v_ice1(ji,jj) = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj) & 697 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & 698 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 699 700 u_ice2(ji,jj) = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1) & 701 & +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) & 702 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 703 ENDIF ! zdummy 586 IF ( vt_i(ji,jj) <= zvmin ) THEN 587 v_ice1(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji, jj-1) ) * e1t(ji+1,jj) & 588 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) & 589 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 590 591 u_ice2(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) & 592 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) & 593 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 594 ENDIF 704 595 END DO 705 596 END DO 706 597 707 CALL lbc_lnk( u_ice2(:,:), 'V', -1. ) 708 CALL lbc_lnk( v_ice1(:,:), 'U', -1. ) 598 CALL lbc_lnk_multi( u_ice2(:,:), 'V', -1., v_ice1(:,:), 'U', -1. ) 709 599 710 600 ! Recompute delta, shear and div, inputs for mechanical redistribution 711 !CDIR NOVERRCHK712 601 DO jj = k_j1+1, k_jpj-1 713 !CDIR NOVERRCHK 714 DO ji = fs_2, jpim1 !RB bug no vect opt due to tmi 715 !- zdd(:,:), zdt(:,:): divergence and tension at centre 602 DO ji = fs_2, jpim1 !RB bug no vect opt due to zmask 603 !- divu_i(:,:), zdt(:,:): divergence and tension at centre 716 604 !- zds(:,:): shear on northeast corner of grid cells 717 zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) ) 718 !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 ) 719 zdummy = vt_i(ji,jj) 720 IF ( zdummy .LE. hminrhg ) THEN 721 722 zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) & 723 & -e2u(ji-1,jj)*u_ice(ji-1,jj) & 724 & +e1v(ji,jj)*v_ice(ji,jj) & 725 & -e1v(ji,jj-1)*v_ice(ji,jj-1) & 726 & ) & 727 & / area(ji,jj) 728 729 zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj) & 730 & -u_ice(ji-1,jj)/e2u(ji-1,jj) & 731 & )*e2t(ji,jj)*e2t(ji,jj) & 732 & -( v_ice(ji,jj)/e1v(ji,jj) & 733 & -v_ice(ji,jj-1)/e1v(ji,jj-1) & 734 & )*e1t(ji,jj)*e1t(ji,jj) & 735 & ) & 736 & / area(ji,jj) 605 IF ( vt_i(ji,jj) <= zvmin ) THEN 606 607 divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj ) * u_ice(ji-1,jj ) & 608 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji ,jj-1) * v_ice(ji ,jj-1) & 609 & ) * r1_e1e2t(ji,jj) 610 611 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) & 612 & -( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 613 & ) * r1_e1e2t(ji,jj) 737 614 ! 738 615 ! SB modif because ocean has no slip boundary condition 739 zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1) & 740 & - u_ice(ji,jj) / e1u(ji,jj) ) & 741 & * e1f(ji,jj) * e1f(ji,jj) & 742 & + ( v_ice(ji+1,jj) / e2v(ji+1,jj) & 743 & - v_ice(ji,jj) / e2v(ji,jj) ) & 744 & * e2f(ji,jj) * e2f(ji,jj) ) & 745 & / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) & 746 & * tmi(ji,jj) * tmi(ji,jj+1) & 747 & * tmi(ji+1,jj) * tmi(ji+1,jj+1) 748 749 zdst(ji,jj) = ( e2u( ji , jj ) * v_ice1(ji ,jj ) & 750 & - e2u( ji-1, jj ) * v_ice1(ji-1,jj ) & 751 & + e1v( ji , jj ) * u_ice2(ji ,jj ) & 752 & - e1v( ji , jj-1 ) * u_ice2(ji ,jj-1) ) / area(ji,jj) 753 754 ! deltat(ji,jj) = SQRT( zdd(ji,jj)*zdd(ji,jj) & 755 ! & + ( zdt(ji,jj)*zdt(ji,jj) + zdst(ji,jj)*zdst(ji,jj) ) * usecc2 & 756 ! & ) + creepl 757 ! MV rewriting 758 delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst(ji,jj)*zdst(ji,jj) ) * usecc2 ) 759 deltat(ji,jj) = delta + creepl 760 ! END MV 616 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) & 617 & +( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 618 & ) * r1_e1e2f(ji,jj) * ( 2.0 - fmask(ji,jj,1) ) & 619 & * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1) 620 621 zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj ) * v_ice1(ji-1,jj ) & 622 & + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji ,jj-1) * u_ice2(ji ,jj-1) ) * r1_e1e2t(ji,jj) 623 624 delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 ) 625 delta_i(ji,jj) = delta + rn_creepl 761 626 762 ENDIF ! zdummy 763 764 END DO !jj 765 END DO !ji 627 ENDIF 628 END DO 629 END DO 766 630 ! 767 631 !------------------------------------------------------------------------------! 768 632 ! 5) Store stress tensor and its invariants 769 633 !------------------------------------------------------------------------------! 770 !771 634 ! * Invariants of the stress tensor are required for limitd_me 772 635 ! (accelerates convergence and improves stability) 773 636 DO jj = k_j1+1, k_jpj-1 774 637 DO ji = fs_2, fs_jpim1 775 divu_i (ji,jj) = zdd (ji,jj) 776 delta_i(ji,jj) = deltat(ji,jj) 777 ! begin TECLIM change 778 zdst(ji,jj)= ( e2u( ji , jj ) * v_ice1(ji,jj) & 779 & - e2u( ji-1, jj ) * v_ice1(ji-1,jj) & 780 & + e1v( ji , jj ) * u_ice2(ji,jj) & 781 & - e1v( ji , jj-1 ) * u_ice2(ji,jj-1) ) / area(ji,jj) 782 shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst(ji,jj) * zdst(ji,jj) ) 783 ! end TECLIM change 638 zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u( ji-1, jj ) * v_ice1(ji-1,jj) & 639 & + e1v(ji,jj) * u_ice2(ji,jj) - e1v( ji , jj-1 ) * u_ice2(ji,jj-1) ) * r1_e1e2t(ji,jj) 640 shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst ) 784 641 END DO 785 642 END DO 786 643 787 644 ! Lateral boundary condition 788 CALL lbc_lnk( divu_i (:,:), 'T', 1. ) 789 CALL lbc_lnk( delta_i(:,:), 'T', 1. ) 790 ! CALL lbc_lnk( shear_i(:,:), 'F', 1. ) 791 CALL lbc_lnk( shear_i(:,:), 'T', 1. ) 645 CALL lbc_lnk_multi( divu_i (:,:), 'T', 1., delta_i(:,:), 'T', 1., shear_i(:,:), 'T', 1. ) 792 646 793 647 ! * Store the stress tensor for the next time step … … 820 674 DO jj = k_j1+1, k_jpj-1 821 675 DO ji = 2, jpim1 822 IF (zpresh(ji,jj) .GT.1.0) THEN676 IF (zpresh(ji,jj) > 1.0) THEN 823 677 sigma1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) 824 678 sigma2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) … … 834 688 ! 835 689 CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 836 CALL wrk_dealloc( jpi,jpj, zc1 , u_oce1, u_oce2, u_ice2, zusw , v_oce1 , v_oce2, v_ice1)837 CALL wrk_dealloc( jpi,jpj, zf1 , deltat, zu_ice, zf2 , deltac, zv_ice , zdd , zdt , zds , zdst)838 CALL wrk_dealloc( jpi,jpj, zd d , zdt , zds , zs1 , zs2 , zs12 , zresr , zpice )690 CALL wrk_dealloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask ) 691 CALL wrk_dealloc( jpi,jpj, zf1 , zu_ice, zf2 , zv_ice , zdt , zds ) 692 CALL wrk_dealloc( jpi,jpj, zdt , zds , zs1 , zs2 , zs12 , zresr , zpice ) 839 693 840 694 END SUBROUTINE lim_rhg
Note: See TracChangeset
for help on using the changeset viewer.
