Changeset 834 for trunk/NEMO/LIM_SRC_3/limrhg.F90
- Timestamp:
- 2008-03-07T18:11:35+01:00 (16 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/LIM_SRC_3/limrhg.F90 (modified) (34 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/LIM_SRC_3/limrhg.F90
r825 r834 2 2 !!====================================================================== 3 3 !! *** MODULE limrhg *** 4 !! Ice rheology : performssea ice rheology4 !! Ice rheology : sea ice rheology 5 5 !!====================================================================== 6 6 #if defined key_lim3 7 7 !!---------------------------------------------------------------------- 8 !! 'key_lim3' LIMsea-ice model8 !! 'key_lim3' LIM3 sea-ice model 9 9 !!---------------------------------------------------------------------- 10 10 !! lim_rhg : computes ice velocities … … 36 36 rone = 1.e0 37 37 !!---------------------------------------------------------------------- 38 !! LIM 2.0, UCL-LOCEAN-IPSL (2005)38 !! LIM 3.0, UCL-LOCEAN-IPSL (2008) 39 39 !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limrhg.F90,v 1.5 2005/03/27 18:34:42 opalod Exp $ 40 40 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt … … 44 44 45 45 SUBROUTINE lim_rhg( k_j1, k_jpj ) 46 46 47 !!------------------------------------------------------------------- 47 !! *** SUBROUTINR lim_rhg *** 48 !! 49 !! ** purpose : determines sea ice drift from wind stress, ice-ocean 48 !! *** SUBROUTINE lim_rhg *** 49 !! EVP-C-grid 50 !! 51 !! ** purpose : determines sea ice drift from wind stress, ice-ocean 50 52 !! stress and sea-surface slope. Ice-ice interaction is described by 51 !! a non-linear elasto-viscous-plastic law including shear strength 52 !! and a bulk rheology (Hunke and Dukowicz, 2002). 53 !! 54 !! Grid B: 55 !! the routine calculates 4 estimates of the divergence, tension, 56 !! and shear on each corner of a given scalar grid box. Likewise, 57 !! four estimates of the components of the stress tensor are 58 !! calculated on each corner. The internal forces on a corner are 59 !! calculated as averages of the four stress tensor contributions. 60 !! 61 !! ** Action : - compute u_ice, v_ice the sea-ice velocity 62 !! 63 !! ** NOTE : for the ice dynamics, the labeling convention is 64 !! that the velocity point (i,j) is located on the southwest corner 65 !! of a scalar (i,j) grid box. This choice is somewhat unfortunate, 66 !! as most models locate the velocity point (i,j) on the northeast 67 !! corner... 68 !! 69 !! ** MORE IMPORTANT NOTES : related to the note above, it is crucial 70 !! to make sure that all variables and coefficients are located on 71 !! right points of the grid. Verify everywhere! 53 !! a non-linear elasto-viscous-plastic (EVP) law including shear 54 !! strength and a bulk rheology (Hunke and Dukowicz, 2002). 55 !! 56 !! The points in the C-grid look like this, dear reader 57 !! 58 !! (ji,jj) 59 !! | 60 !! | 61 !! (ji-1,jj) | (ji,jj) 62 !! --------- 63 !! | | 64 !! | (ji,jj) |------(ji,jj) 65 !! | | 66 !! --------- 67 !! (ji-1,jj-1) (ji,jj-1) 68 !! 69 !! ** Inputs : - wind forcing (stress), oceanic currents 70 !! ice total volume (vt_i) per unit area 71 !! snow total volume (vt_s) per unit area 72 !! 73 !! ** Action : - compute u_ice, v_ice : the components of the 74 !! sea-ice velocity vector 75 !! - compute delta_i, shear_i, divu_i, which are inputs 76 !! of the ice thickness distribution 77 !! 78 !! ** Steps : 1) Compute ice snow mass, ice strength 79 !! 2) Compute wind, oceanic stresses, mass terms and 80 !! coriolis terms of the momentum equation 81 !! 3) Solve the momentum equation (iterative procedure) 82 !! 4) Prevent high velocities if the ice is thin 83 !! 5) Recompute invariants of the strain rate tensor 84 !! which are inputs of the ITD, store stress 85 !! for the next time step 86 !! 6) Control prints of residual (convergence) 87 !! and charge ellipse. 88 !! The user should make sure that the parameters 89 !! nevp, telast and creepl maintain stress state 90 !! on the charge ellipse for plastic flow 91 !! e.g. in the Canadian Archipelago 92 !! 93 !! ** References : Hunke and Dukowicz, JPO97 94 !! Bouillon et al., 08, in prep (update this when 95 !! published) 96 !! Vancoppenolle et al., OM08 72 97 !! 73 98 !! History : 74 !! 1.0 ! 07-03 (M.A. Morales Maqueda, S. Bouillon , M. Vancoppenolle)75 !! EVP, C grid,LIM399 !! 1.0 ! 07-03 (M.A. Morales Maqueda, S. Bouillon) 100 !! 2.0 ! 08-03 M. Vancoppenolle : LIM3 76 101 !! 77 102 !!------------------------------------------------------------------- … … 79 104 ! 80 105 INTEGER, INTENT(in) :: & 81 k_j1 , &!: southern j-index for ice computation82 k_jpj !: northern j-index for ice computation106 k_j1 , & !: southern j-index for ice computation 107 k_jpj !: northern j-index for ice computation 83 108 84 109 ! * Local variables 85 INTEGER :: ji, jj , jl!: dummy loop indices110 INTEGER :: ji, jj !: dummy loop indices 86 111 87 112 INTEGER :: & 88 iim1, ijm1, iip1 , ijp1 , & ! temporary integers 89 iter, jter ! " " 113 iter, jter !: temporary integers 90 114 91 115 CHARACTER (len=50) :: charout 92 116 93 117 REAL(wp) :: & 94 zt11 , zt12 , zt21 , zt22 , & ! temporary scalars 95 ztagnx, ztagny , delta ! " " 118 zt11, zt12, zt21, zt22, & !: temporary scalars 119 ztagnx, ztagny, & !: wind stress on U/V points 120 delta ! 96 121 97 122 REAL(wp) :: & 98 za, zstms, zsang, zmask 99 100 REAL(wp),DIMENSION(jpi,jpj) :: & 101 zpresh, zpreshc, zfrld1, zfrld2, zmass1, zmass2, zcorl1, zcorl2, & 102 za1ct, za2ct, zc1, zusw , & 103 u_oce1, v_oce1, u_oce2, v_oce2, u_ice2, v_ice1 123 za, & !: 124 zstms, & !: temporary scalar for ice strength 125 zsang, & !: temporary scalar for coriolis term 126 zmask !: mask for the computation of ice mass 127 128 REAL(wp),DIMENSION(jpi,jpj) :: & 129 zpresh , & !: temporary array for ice strength 130 zpreshc , & !: Ice strength on grid cell corners (zpreshc) 131 zfrld1, zfrld2, & !: lead fraction on U/V points 132 zmass1, zmass2, & !: ice/snow mass on U/V points 133 zcorl1, zcorl2, & !: coriolis parameter on U/V points 134 za1ct, za2ct , & !: temporary arrays 135 zc1 , & !: ice mass 136 zusw , & !: temporary weight for the computation 137 !: of ice strength 138 u_oce1, v_oce1, & !: ocean u/v component on U points 139 u_oce2, v_oce2, & !: ocean u/v component on V points 140 u_ice2, & !: ice u component on V point 141 v_ice1 !: ice v component on U point 104 142 105 143 REAL(wp) :: & 106 dtevp,dtotel,ecc2,z0,zr,zcca,zccb,denr, & 107 zu_ice2,zv_ice1, zddc, zdtc, zdst, zdsshx, zdsshy, & 108 sigma1, sigma2 109 110 REAL(wp),DIMENSION(jpi,jpj) :: & 111 zf1, zf2 112 113 REAL(wp),DIMENSION(jpi,jpj) :: & 114 zdd, & !: Divergence [ 115 zdt, & !: 116 zds, & !: 117 deltat, & 118 deltac, & 119 zs1, & 120 zs2, & 121 zs12 122 123 REAL :: & 124 zumax !: Maximal tolerated ice velocity 125 126 REAL(wp) :: & 127 zresm !: Maximal error on ice velocity 128 129 REAL(wp),DIMENSION(jpi,jpj) :: & 130 zu_ice , & !: Ice velocity on previous time step 131 zv_ice , & 132 zresr !: Local error on velocity 144 dtevp, & !: time step for subcycling 145 dtotel, & !: 146 ecc2, & !: square of yield ellipse eccenticity 147 z0, & !: temporary scalar 148 zr, & !: temporary scalar 149 zcca, zccb, & !: temporary scalars 150 zu_ice2, & !: 151 zv_ice1, & !: 152 zddc, zdtc, & !: temporary array for delta on corners 153 zdst, & !: temporary array for delta on centre 154 zdsshx, zdsshy, & !: term for the gradient of ocean surface 155 sigma1, sigma2 !: internal ice stress 156 157 REAL(wp),DIMENSION(jpi,jpj) :: & 158 zf1, zf2 !: arrays for internal stresses 159 160 REAL(wp),DIMENSION(jpi,jpj) :: & 161 zdd, zdt, & !: Divergence and tension at centre of grid cells 162 zds, & !: Shear on northeast corner of grid cells 163 deltat, & !: Delta at centre of grid cells 164 deltac, & !: Delta on corners 165 zs1, zs2, & !: Diagonal stress tensor components zs1 and zs2 166 zs12 !: Non-diagonal stress tensor component zs12 133 167 134 168 REAL(wp) :: & 135 zindb , & !: ice (1) or not (0) 136 zdummy !: ice computation 169 zresm , & !: Maximal error on ice velocity 170 zindb , & !: ice (1) or not (0) 171 zdummy !: dummy argument 172 173 REAL(wp),DIMENSION(jpi,jpj) :: & 174 zu_ice , & !: Ice velocity on previous time step 175 zv_ice , & 176 zresr !: Local error on velocity 137 177 138 178 INTEGER :: & 139 179 stress_mean_swi = 1 140 141 !!---------------------------------------------------------------------- 142 143 ! 144 !------------------------------------------------------------------------------! 145 ! 0) Ice-Snow mass (zc1), ice strength (zpresh) ! 146 !------------------------------------------------------------------------------! 147 ! 148 ! Maximal tolerated velocity 149 ! zumax = 1.0 150 180 ! 181 !------------------------------------------------------------------------------! 182 ! 1) Ice-Snow mass (zc1), ice strength (zpresh) ! 183 !------------------------------------------------------------------------------! 184 ! 151 185 ! Put every vector to 0 152 186 zpresh(:,:) = 0.0 ; zpreshc(:,:) = 0.0 ; zfrld1(:,:) = 0.0 … … 154 188 za1ct(:,:) = 0.0 ; za2ct(:,:) = 0.0 155 189 zc1(:,:) = 0.0 ; zusw(:,:) = 0.0 156 157 190 u_oce1(:,:) = 0.0 ; v_oce1(:,:) = 0.0 ; u_oce2(:,:) = 0.0 158 191 v_oce2(:,:) = 0.0 ; u_ice2(:,:) = 0.0 ; v_ice1(:,:) = 0.0 159 160 zf1(:,:) = 0.0 ; zf2(:,:) = 0.0 161 192 zf1(:,:) = 0.0 ; zf2(:,:) = 0.0 162 193 zdd(:,:) = 0.0 ; zdt(:,:) = 0.0 ; zds(:,:) = 0.0 163 194 deltat(:,:) = 0.0 ; deltac(:,:) = 0.0 164 165 ! 166 !------------------------------------------------------------------------------! 167 ! 1) Ice-Snow mass (zc1), ice strength (zpresh) ! 168 !------------------------------------------------------------------------------! 169 ! 170 ! compute ice strength 195 zpresh(:,:) = 0.0 196 197 ! Ice strength on T-points 171 198 CALL lim_itd_me_icestrength(ridge_scheme_swi) 172 199 173 zpresh(:,:) = 0.0 174 200 ! Ice mass and temp variables 175 201 DO jj = k_j1 , k_jpj-1 176 202 DO ji = 1 , jpi … … 178 204 zpresh(ji,jj) = tms(ji,jj) * strength(ji,jj) / 2. 179 205 ! tmi = 1 where ther is ice or on land 206 ! Claude sees a bug, next line : tmi(ji,jj) 180 207 tmi = 1.0 - ( 1.0 - MAX( 0.0 , SIGN ( 1.0 , vt_i(ji,jj) - & 181 208 epsd ) ) ) * tms(ji,jj) … … 186 213 CALL lbc_lnk( zpresh(:,:), 'T', 1. ) 187 214 188 ! Ice strength (zpreshc) on grid cell corners (needed for 189 ! calculation of shear stress 215 ! Claude sees a bug 216 ! CALL lbc_lnk( tmi(:,:), 'T', 1. ) 217 218 ! Ice strength on grid cell corners (zpreshc) 219 ! needed for calculation of shear stress 190 220 DO jj = k_j1+1, k_jpj-1 191 221 DO ji = 2, jpim1 192 zstms = tms(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + tms(ji,jj+1) * wght(ji+1,jj+1,1,2) & 193 & + tms(ji+1,jj ) * wght(ji+1,jj+1,2,1) + tms(ji,jj ) * wght(ji+1,jj+1,1,1) 222 zstms = tms(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + & 223 & tms(ji,jj+1) * wght(ji+1,jj+1,1,2) + & 224 & tms(ji+1,jj) * wght(ji+1,jj+1,2,1) + & 225 & tms(ji,jj) * wght(ji+1,jj+1,1,1) 194 226 zusw(ji,jj) = 1.0 / MAX( zstms, epsd ) 195 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) & 196 & + zpresh(ji+1,jj ) * wght(ji+1,jj+1,2,1) + zpresh(ji,jj ) * wght(ji+1,jj+1,1,1) ) & 197 & * zusw(ji,jj) 227 zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + & 228 & zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) + & 229 & zpresh(ji+1,jj) * wght(ji+1,jj+1,2,1) + & 230 & zpresh(ji,jj) * wght(ji+1,jj+1,1,1) & 231 & ) * zusw(ji,jj) 198 232 END DO 199 233 END DO … … 233 267 234 268 ! Mass, coriolis coeff. and currents 235 zmass1(ji,jj) = (zt12*zc1(ji,jj)+zt11*zc1(ji+1,jj))/(zt11+zt12+epsd) 236 zmass2(ji,jj) = (zt22*zc1(ji,jj)+zt21*zc1(ji,jj+1))/(zt21+zt22+epsd) 237 238 zcorl1(ji,jj) = zmass1(ji,jj)*(e1t(ji+1,jj)*fcor(ji,jj)+e1t(ji,jj)*fcor(ji+1,jj))/(e1t(ji,jj)+e1t(ji+1,jj)+epsd) 239 zcorl2(ji,jj) = zmass2(ji,jj)*(e2t(ji,jj+1)*fcor(ji,jj)+e2t(ji,jj)*fcor(ji,jj+1))/(e2t(ji,jj+1)+e2t(ji,jj)+epsd) 240 ! 241 ! Reminder: since this is a C grid, the location of ocean currents 242 ! calculated by OPA should coincide with ice model vector points: 243 ! no need for interpolation once the definition of u_io and v_io 244 ! has been changed in module "icestp". Arrays u_oce1 and v_oce2 could 245 ! then be replaced by u_oce and v_oce, respectively. 269 zmass1(ji,jj) = ( zt12*zc1(ji,jj) + zt11*zc1(ji+1,jj) ) / (zt11+zt12+epsd) 270 zmass2(ji,jj) = ( zt22*zc1(ji,jj) + zt21*zc1(ji,jj+1) ) / (zt21+zt22+epsd) 271 zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj)*fcor(ji,jj) + & 272 e1t(ji,jj)*fcor(ji+1,jj) ) & 273 / (e1t(ji,jj) + e1t(ji+1,jj) + epsd ) 274 zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1)*fcor(ji,jj) + & 275 e2t(ji,jj)*fcor(ji,jj+1) ) & 276 / ( e2t(ji,jj+1) + e2t(ji,jj) + epsd ) 246 277 ! 247 278 u_oce1(ji,jj) = u_io(ji,jj) 248 249 ! SB modif because ocean has no slip boundary condition 279 v_oce2(ji,jj) = v_io(ji,jj) 280 281 ! Ocean has no slip boundary condition 250 282 v_oce1(ji,jj) = 0.5*( (v_io(ji,jj)+v_io(ji,jj-1))*e1t(ji+1,jj) & 251 283 & +(v_io(ji+1,jj)+v_io(ji+1,jj-1))*e1t(ji,jj)) & 252 284 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 253 285 254 u_oce2(ji,jj) = 0.5*( (u_io(ji,jj)+u_io(ji-1,jj))*e2t(ji,jj+1)&286 u_oce2(ji,jj) = 0.5*((u_io(ji,jj)+u_io(ji-1,jj))*e2t(ji,jj+1) & 255 287 & +(u_io(ji,jj+1)+u_io(ji-1,jj+1))*e2t(ji,jj)) & 256 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 257 258 v_oce2(ji,jj) = v_io(ji,jj) 288 & / (e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 259 289 260 290 ! Wind stress. … … 262 292 ztagny = ( 1. - zfrld2(ji,jj) ) * gtauy(ji,jj) 263 293 264 ! Computation of the velocity field taking into account the ice -ice interaction.294 ! Computation of the velocity field taking into account the ice internal interaction. 265 295 ! Terms that are independent of the velocity field. 266 ! Reminder: does zcorl*(-v_oce,u_oce) represent the surface pressure gradient? Why is it still267 ! calculated in this way? An ocean model with a free surface should provide the gradient268 ! of surface elevation directly...269 296 270 297 ! SB On utilise maintenant le gradient de la pente de l'ocean … … 272 299 ! zdsshx = (ssh_io(ji+1,jj) - ssh_io(ji,jj))/e1u(ji,jj) 273 300 ! zdsshy = (ssh_io(ji,jj+1) - ssh_io(ji,jj))/e2v(ji,jj) 301 274 302 zdsshx = 0.0 275 303 zdsshy = 0.0 … … 283 311 ! 284 312 !------------------------------------------------------------------------------! 285 ! 3) Solution of the momentum equation 313 ! 3) Solution of the momentum equation, iterative procedure 286 314 !------------------------------------------------------------------------------! 287 315 ! … … 319 347 END DO 320 348 321 322 349 DO jj = k_j1+1, k_jpj-1 350 DO ji = 2, jpim1 323 351 324 352 ! … … 327 355 !- zds(:,:): shear on northeast corner of grid cells 328 356 ! 329 !- IMPORTANT REMINDER: note that, the way these terms are coded, there are many repeated 330 ! calculations. Speed could be improved by regrouping terms. For 357 !- IMPORTANT REMINDER: Dear Gurvan, note that, the way these terms are coded, 358 ! there are many repeated calculations. 359 ! Speed could be improved by regrouping terms. For 331 360 ! the moment, however, the stress is on clarity of coding to avoid 332 ! bugs (mamm). 333 !- ALSO: arrays zdd, zdt, zds and delta could be removed in the future to minimise memory demand. 361 ! bugs (Martin, for Miguel). 362 ! 363 !- ALSO: arrays zdd, zdt, zds and delta could 364 ! be removed in the future to minimise memory demand. 365 ! 334 366 !- MORE NOTES: Note that we are calculating deformation rates and stresses on the corners of 335 367 ! grid cells, exactly as in the B grid case. For simplicity, the indexation on … … 337 369 ! 338 370 ! 339 ! (ji,jj)340 ! |341 ! |342 ! (ji-1,jj) | (ji,jj)343 ! ---------344 ! | |345 ! | (ji,jj) |------(ji,jj)346 ! | |347 ! ---------348 !(ji-1,jj-1) (ji,jj-1)349 !350 351 371 zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) & 352 372 & -e2u(ji-1,jj)*u_ice(ji-1,jj) & … … 366 386 367 387 ! 368 ! SB modif because ocean has no slip boundary condition369 388 zds(ji,jj) = ( ( u_ice(ji,jj+1)/e1u(ji,jj+1) & 370 389 & -u_ice(ji,jj)/e1u(ji,jj) & … … 379 398 380 399 381 ! SB modif because need to compute zddc and zdtc correctly382 400 v_ice1(ji,jj) = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj) & 383 401 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & … … 406 424 & / area(ji,jj) 407 425 408 ! Old lines409 ! deltat(ji,jj) = SQRT( zdd(ji,jj)*zdd(ji,jj) + &410 ! & ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 &411 ! & ) + creepl412 ! New lines413 ! Regularisation of dP/dx414 426 delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + & 415 427 & ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 ) 416 428 deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 + & 417 429 (zdt(ji,jj)**2 + zdst**2)*usecc2), creepl ) 418 ! End of new lines419 430 420 431 !-Calculate stress tensor components zs1 and zs2 421 432 !-at centre of grid cells (see section 3.5 of CICE user's guide). 422 ! Old lines423 ! zs1(ji,jj) = ( zs1(ji,jj) &424 ! & - dtotel*((1.0-alphaevp)*zs1(ji,jj)+(1.0-zdd(ji,jj)/deltat(ji,jj))*zpresh(ji,jj))) &425 ! & /(1.0+alphaevp*dtotel)426 ! New lines427 433 zs1(ji,jj) = ( zs1(ji,jj) & 428 434 & - dtotel*( ( 1.0 - alphaevp) * zs1(ji,jj) + & 429 & ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) * zpresh(ji,jj) ) ) & 435 & ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) & 436 * zpresh(ji,jj) ) ) & 430 437 & / ( 1.0 + alphaevp * dtotel ) 431 ! End of new lines 432 433 zs2(ji,jj) = ( zs2(ji,jj)&434 & - dtotel*((1.0-alphaevp)*ecc2*zs2(ji,jj)-zdt(ji,jj)/deltat(ji,jj)*zpresh(ji,jj)))&435 & / (1.0+alphaevp*ecc2*dtotel)438 439 zs2(ji,jj) = ( zs2(ji,jj) & 440 & - dtotel*((1.0-alphaevp)*ecc2*zs2(ji,jj) - & 441 zdt(ji,jj)/deltat(ji,jj)*zpresh(ji,jj)) ) & 442 & / ( 1.0 + alphaevp*ecc2*dtotel ) 436 443 437 444 END DO … … 443 450 DO jj = k_j1+1, k_jpj-1 444 451 DO ji = 2, jpim1 445 446 452 !- Calculate Delta on corners 447 448 453 zddc = ( ( v_ice1(ji,jj+1)/e1u(ji,jj+1) & 449 454 & -v_ice1(ji,jj)/e1u(ji,jj) & … … 455 460 & / ( e1f(ji,jj) * e2f(ji,jj) ) 456 461 457 458 462 zdtc = (-( v_ice1(ji,jj+1)/e1u(ji,jj+1) & 459 463 & -v_ice1(ji,jj)/e1u(ji,jj) & … … 465 469 & / ( e1f(ji,jj) * e2f(ji,jj) ) 466 470 467 deltac(ji,jj) = sqrt(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl471 deltac(ji,jj) = SQRT(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl 468 472 469 473 !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide). 470 471 zs12(ji,jj) = ( zs12(ji,jj)&472 & -dtotel*((1.0-alphaevp)*ecc2*zs12(ji,jj)-zds(ji,jj)/(2.0*deltac(ji,jj))*zpreshc(ji,jj)))&473 & / (1.0+alphaevp*ecc2*dtotel)474 475 END DO 476 END DO 474 zs12(ji,jj) = ( zs12(ji,jj) & 475 & - dtotel*( (1.0-alphaevp)*ecc2*zs12(ji,jj) - zds(ji,jj) / & 476 & ( 2.0*deltac(ji,jj) ) * zpreshc(ji,jj))) & 477 & / ( 1.0 + alphaevp*ecc2*dtotel ) 478 479 END DO ! ji 480 END DO ! jj 477 481 478 482 CALL lbc_lnk( zs12(:,:), 'F', 1. ) … … 481 485 DO jj = k_j1+1, k_jpj-1 482 486 DO ji = 2, jpim1 483 484 ! 485 !- contribution of zs1, zs2 and zs12 to zf1 486 ! 487 ! (ji,jj) 488 ! | 489 ! | 490 ! (ji-1,jj) | (ji,jj) 491 ! --------- 492 ! | | 493 ! | (ji,jj) |------(ji,jj) 494 ! | | 495 ! --------- 496 !(ji-1,jj-1) (ji,jj-1) 497 ! 498 499 zf1(ji,jj) = 0.5*( (zs1(ji+1,jj)-zs1(ji,jj))*e2u(ji,jj) & 500 & +(zs2(ji+1,jj)*e2t(ji+1,jj)**2-zs2(ji,jj)*e2t(ji,jj)**2)/e2u(ji,jj) & 501 & +2.0*(zs12(ji,jj)*e1f(ji,jj)**2-zs12(ji,jj-1)*e1f(ji,jj-1)**2)/e1u(ji,jj) & 502 & ) & 503 & /(e1u(ji,jj)*e2u(ji,jj)) 504 505 ! 506 !contribution of zs1, zs2 and zs12 to zf2 507 ! 508 ! (ji,jj) 509 ! | 510 ! | 511 ! (ji-1,jj) | (ji,jj) 512 ! --------- 513 ! | | 514 ! | (ji,jj) |------(ji,jj) 515 ! | | 516 ! --------- 517 !(ji-1,jj-1) (ji,jj-1) 518 ! 519 520 zf2(ji,jj) = 0.5*( (zs1(ji,jj+1)-zs1(ji,jj))*e1v(ji,jj) & 521 & -(zs2(ji,jj+1)*e1t(ji,jj+1)**2-zs2(ji,jj)*e1t(ji,jj)**2)/e1v(ji,jj) & 522 & +2.0*(zs12(ji,jj)*e2f(ji,jj)**2-zs12(ji-1,jj)*e2f(ji-1,jj)**2)/e2v(ji,jj) & 523 & ) & 524 & /(e1v(ji,jj)*e2v(ji,jj)) 525 487 !- contribution of zs1, zs2 and zs12 to zf1 488 zf1(ji,jj) = 0.5*( (zs1(ji+1,jj)-zs1(ji,jj))*e2u(ji,jj) & 489 & +(zs2(ji+1,jj)*e2t(ji+1,jj)**2-zs2(ji,jj)*e2t(ji,jj)**2)/e2u(ji,jj) & 490 & +2.0*(zs12(ji,jj)*e1f(ji,jj)**2-zs12(ji,jj-1)*e1f(ji,jj-1)**2)/e1u(ji,jj) & 491 & ) / ( e1u(ji,jj)*e2u(ji,jj) ) 492 ! contribution of zs1, zs2 and zs12 to zf2 493 zf2(ji,jj) = 0.5*( (zs1(ji,jj+1)-zs1(ji,jj))*e1v(ji,jj) & 494 & -(zs2(ji,jj+1)*e1t(ji,jj+1)**2 - zs2(ji,jj)*e1t(ji,jj)**2)/e1v(ji,jj) & 495 & + 2.0*(zs12(ji,jj)*e2f(ji,jj)**2 - & 496 zs12(ji-1,jj)*e2f(ji-1,jj)**2)/e2v(ji,jj) & 497 & ) / ( e1v(ji,jj)*e2v(ji,jj) ) 526 498 END DO 527 499 END DO … … 531 503 ! Both the Coriolis term and the ice-ocean drag are solved semi-implicitly. 532 504 ! 533 IF ( mod(iter,2).eq.0) THEN505 IF (MOD(iter,2).eq.0) THEN 534 506 535 507 DO jj = k_j1+1, k_jpj-1 536 508 DO ji = 2, jpim1 537 ! 538 ! (ji,jj) 539 ! | 540 ! | 541 ! (ji-1,jj) | (ji,jj) 542 ! --------- 543 ! | | 544 ! | (ji,jj) |------(ji,jj) 545 ! | | 546 ! --------- 547 !(ji-1,jj-1) (ji,jj-1) 548 ! 549 zmask = (1.0-max(rzero,sign(rone,-zmass1(ji,jj))))*tmu(ji,jj) 509 zmask = (1.0-MAX(rzero,SIGN(rone,-zmass1(ji,jj))))*tmu(ji,jj) 550 510 zsang = SIGN ( 1.0 , fcor(ji,jj) ) * sangvg 551 511 z0 = zmass1(ji,jj)/dtevp … … 555 515 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & 556 516 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 557 za = rhoco*sqrt((u_ice(ji,jj)-u_oce1(ji,jj))**2+(zv_ice1-v_oce1(ji,jj))**2)*(1.0-zfrld1(ji,jj)) 558 zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za*(cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj)) 517 za = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + & 518 (zv_ice1-v_oce1(ji,jj))**2) * (1.0-zfrld1(ji,jj)) 519 zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + & 520 za*(cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj)) 559 521 zcca = z0+za*cangvg 560 522 zccb = zcorl1(ji,jj)+za*zsang 561 523 u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask 562 ! u_ice(ji,jj) = MAX( -1.0 , MIN( u_ice(ji,jj), 1.0 ) )563 524 564 525 END DO … … 569 530 DO jj = k_j1+1, k_jpj-1 570 531 DO ji = 2, jpim1 571 ! 572 ! (ji,jj) 573 ! | 574 ! | 575 ! (ji-1,jj) | (ji,jj) 576 ! --------- 577 ! | | 578 ! | (ji,jj) |------(ji,jj) 579 ! | | 580 ! --------- 581 !(ji-1,jj-1) (ji,jj-1) 582 ! 532 583 533 zmask = (1.0-max(rzero,sign(rone,-zmass2(ji,jj))))*tmv(ji,jj) 584 534 zsang = sign(1.0,fcor(ji,jj))*sangvg … … 588 538 & + (u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) & 589 539 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 590 za = rhoco*sqrt((zu_ice2-u_oce2(ji,jj))**2+(v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 591 zr = z0*v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za*(cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj)) 540 za = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + & 541 (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 542 zr = z0*v_ice(ji,jj) + zf2(ji,jj) + & 543 za2ct(ji,jj) + za*(cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj)) 592 544 zcca = z0+za*cangvg 593 545 zccb = zcorl2(ji,jj)+za*zsang 594 546 v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask 595 ! v_ice(ji,jj) = MAX( -1.0 , MIN( v_ice(ji,jj), 1.0 ) )596 547 597 548 END DO … … 600 551 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 601 552 602 553 ELSE 603 554 DO jj = k_j1+1, k_jpj-1 604 555 DO ji = 2, jpim1 605 !606 ! (ji,jj)607 ! |608 ! |609 ! (ji-1,jj) | (ji,jj)610 ! ---------611 ! | |612 ! | (ji,jj) |------(ji,jj)613 ! | |614 ! ---------615 !(ji-1,jj-1) (ji,jj-1)616 !617 556 zmask = (1.0-max(rzero,sign(rone,-zmass2(ji,jj))))*tmv(ji,jj) 618 557 zsang = sign(1.0,fcor(ji,jj))*sangvg … … 623 562 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 624 563 625 za = rhoco*sqrt((zu_ice2-u_oce2(ji,jj))**2+(v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 626 zr = z0*v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za*(cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj)) 564 za = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + & 565 (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 566 zr = z0*v_ice(ji,jj) + zf2(ji,jj) + & 567 za2ct(ji,jj) + za*(cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj)) 627 568 zcca = z0+za*cangvg 628 569 zccb = zcorl2(ji,jj)+za*zsang 629 570 v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask 630 ! v_ice(ji,jj) = MAX( -1.0 , MIN( v_ice(ji,jj), 1.0 ) )631 571 632 572 END DO … … 637 577 DO jj = k_j1+1, k_jpj-1 638 578 DO ji = 2, jpim1 639 ! 640 ! (ji,jj) 641 ! | 642 ! | 643 ! (ji-1,jj) | (ji,jj) 644 ! --------- 645 ! | | 646 ! | (ji,jj) |------(ji,jj) 647 ! | | 648 ! --------- 649 ! (ji-1,jj-1) (ji,jj-1) 650 ! 651 652 zmask = (1.0-max(rzero,sign(rone,-zmass1(ji,jj))))*tmu(ji,jj) 653 zsang = sign(1.0,fcor(ji,jj))*sangvg 579 zmask = (1.0-MAX(rzero,SIGN(rone,-zmass1(ji,jj))))*tmu(ji,jj) 580 zsang = SIGN(1.0,fcor(ji,jj))*sangvg 654 581 z0 = zmass1(ji,jj)/dtevp 655 582 ! SB modif because ocean has no slip boundary condition … … 658 585 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 659 586 660 za = rhoco*sqrt((u_ice(ji,jj)-u_oce1(ji,jj))**2+(zv_ice1-v_oce1(ji,jj))**2)*(1.0-zfrld1(ji,jj)) 661 zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za*(cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj)) 587 za = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + & 588 (zv_ice1-v_oce1(ji,jj))**2)*(1.0-zfrld1(ji,jj)) 589 zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + & 590 za*(cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj)) 662 591 zcca = z0+za*cangvg 663 592 zccb = zcorl1(ji,jj)+za*zsang 664 593 u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask 665 ! u_ice(ji,jj) = MAX( -1.0 , MIN( u_ice(ji,jj), 1.0 ) ) 666 667 END DO 668 END DO 594 END DO ! ji 595 END DO ! jj 669 596 670 597 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 671 598 672 ENDIF 673 674 !--- Convergence test. 675 DO jj = k_j1+1 , k_jpj-1 676 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) , & 677 ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 678 END DO 679 zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) ) 680 681 !------------------------------------------- 682 ! Compute dynamical inputs of the itd model 683 !------------------------------------------- 684 ! Mean over all iterations 685 686 ! IF ( stress_mean_swi .EQ. 1 ) THEN 687 ! DO jj = k_j1+1, k_jpj-1 688 ! DO ji = 2, jpim1 689 ! divu_i(ji,jj) = divu_i(ji,jj) + zdd(ji,jj) / nevp 690 ! delta_i(ji,jj) = delta_i(ji,jj) + deltat (ji,jj) / nevp 691 ! shear_i(ji,jj) = shear_i(ji,jj) + zds (ji,jj) / nevp 692 ! END DO 693 ! END DO 694 ! ENDIF 599 ENDIF 600 601 !--- Convergence test. 602 DO jj = k_j1+1 , k_jpj-1 603 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) , & 604 ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 605 END DO 606 zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) ) 695 607 696 608 ! ! ==================== ! … … 698 610 ! ! ==================== ! 699 611 700 !------------------------- 701 ! Prevent high velocities 702 !------------------------- 612 ! 613 !------------------------------------------------------------------------------! 614 ! 4) Prevent ice velocities when the ice is thin 615 !------------------------------------------------------------------------------! 616 ! 703 617 ! If the ice thickness is below 1cm then ice velocity should equal the 704 ! ocean velocity 618 ! ocean velocity, 619 ! This prevents high velocity when ice is thin 705 620 DO jj = k_j1+1, k_jpj-1 706 621 DO ji = 2, jpim1 … … 713 628 END DO 714 629 END DO 715 630 ! 631 !------------------------------------------------------------------------------! 632 ! 5) Compute stress rain invariants and store strain tensor 633 !------------------------------------------------------------------------------! 634 ! 635 ! Compute delta, shear and div, inputs for mechanical redistribution 716 636 DO jj = k_j1+1, k_jpj-1 717 637 DO ji = 2, jpim1 718 ! Recompute stress tensor invariants719 638 !- zdd(:,:), zdt(:,:): divergence and tension at centre 720 ! of grid cells721 639 !- zds(:,:): shear on northeast corner of grid cells 722 723 640 zindb = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) ) 724 641 zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) … … 782 699 END DO !ji 783 700 784 ! dynamical invariants 701 ! MV This should be removed as it is the same as previous lines 702 ! ! dynamical invariants 785 703 ! IF ( stress_mean_swi .EQ. 0 ) THEN 704 ! the following should not be necessary 786 705 DO jj = k_j1+1, k_jpj-1 787 706 DO ji = 2, jpim1 … … 797 716 CALL lbc_lnk( shear_i(:,:), 'F', 1. ) 798 717 799 IF(ln_ctl) THEN 800 WRITE(charout,FMT="('lim_rhg : res =',D23.16, ' iter =',I4)") zresm, jter 801 CALL prt_ctl_info(charout) 802 CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :') 803 ENDIF 804 805 ! Stress tensor 718 719 ! Store the stress tensor for next time step 720 ! this accelerates convergence and improves stability 806 721 IF ( stress_mean_swi .EQ. 1 ) THEN 807 722 DO jj = k_j1+1, k_jpj-1 … … 814 729 ENDIF 815 730 816 ! Ice internal stresses731 ! Lateral boundary condition 817 732 CALL lbc_lnk( stress1_i(:,:), 'T', 1. ) 818 733 CALL lbc_lnk( stress2_i(:,:), 'T', 1. ) 819 734 CALL lbc_lnk( stress12_i(:,:), 'F', 1. ) 820 735 821 !------------------------ 822 ! Write charge ellipse 823 !------------------------ 824 736 ! 737 !------------------------------------------------------------------------------! 738 ! 6) Control prints of residual and charge ellipse 739 !------------------------------------------------------------------------------! 740 ! 741 ! print the residual for convergence 742 IF(ln_ctl) THEN 743 WRITE(charout,FMT="('lim_rhg : res =',D23.16, ' iter =',I4)") zresm, iter 744 CALL prt_ctl_info(charout) 745 CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :') 746 ENDIF 747 748 ! print charge ellipse 749 ! This can be desactivated once the user is sure that the stress state 750 ! lie on the charge ellipse. See Bouillon et al. 08 for more details 825 751 IF(ln_ctl) THEN 826 752 CALL prt_ctl_info('lim_rhg : numit :',ivar1=numit)
Note: See TracChangeset
for help on using the changeset viewer.
