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