Changeset 2528 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
- Timestamp:
- 2010-12-27T18:33:53+01:00 (13 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
r2477 r2528 7 7 !! 3.0 ! 2008-03 (M. Vancoppenolle) LIM3 8 8 !! - ! 2008-11 (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy 9 !! 3.3 ! 2009-05 (G.Garric) addition of the lim2_evp cas 9 10 !!---------------------------------------------------------------------- 10 #if defined key_lim3 11 #if defined key_lim3 || ( defined key_lim2 && ! defined key_lim2_vp ) 11 12 !!---------------------------------------------------------------------- 12 !! 'key_lim3' LIM3 sea-ice model 13 !! 'key_lim3' OR LIM-3 sea-ice model 14 !! 'key_lim2' AND NOT 'key_lim2_vp' VP LIM-2 sea-ice model 13 15 !!---------------------------------------------------------------------- 14 16 !! lim_rhg : computes ice velocities 15 17 !!---------------------------------------------------------------------- 16 !! * Modules used 17 USE phycst 18 USE par_oce 19 USE dom_oce 20 USE dom_ice 21 USE sbc_oce ! Surface boundary condition: ocean fields 22 USE sbc_ice ! Surface boundary condition: ice fields 23 USE ice 24 USE lbclnk 25 USE lib_mpp 26 USE in_out_manager ! I/O manager 27 USE limitd_me 28 USE prtctl ! Print control 29 18 USE phycst ! Physical constant 19 USE par_oce ! Ocean parameters 20 USE dom_oce ! Ocean domain 21 USE sbc_oce ! Surface boundary condition: ocean fields 22 USE sbc_ice ! Surface boundary condition: ice fields 23 USE lbclnk ! Lateral Boundary Condition / MPP link 24 USE lib_mpp ! MPP library 25 USE in_out_manager ! I/O manager 26 USE prtctl ! Print control 27 #if defined key_lim3 28 USE ice ! LIM-3: ice variables 29 USE dom_ice ! LIM-3: ice domain 30 USE limitd_me ! LIM-3: 31 #else 32 USE ice_2 ! LIM2: ice variables 33 USE dom_ice_2 ! LIM2: ice domain 34 #endif 30 35 31 36 IMPLICIT NONE 32 37 PRIVATE 33 38 34 !! * Routine accessibility 35 PUBLIC lim_rhg ! routine called by lim_dyn 36 39 PUBLIC lim_rhg ! routine called by lim_dyn (or lim_dyn_2) 40 41 REAL(wp) :: rzero = 0._wp ! constant values 42 REAL(wp) :: rone = 1._wp ! constant values 43 37 44 !! * Substitutions 38 45 # include "vectopt_loop_substitute.h90" 39 40 !! * Module variables41 REAL(wp) :: & ! constant values42 rzero = 0.e0 , &43 rone = 1.e044 46 !!---------------------------------------------------------------------- 45 !! LIM 3.0, UCL-LOCEAN-IPSL (2008)47 !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010) 46 48 !! $Id$ 47 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)49 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 48 50 !!---------------------------------------------------------------------- 49 50 51 CONTAINS 51 52 52 53 SUBROUTINE lim_rhg( k_j1, k_jpj ) 53 54 54 !!------------------------------------------------------------------- 55 55 !! *** SUBROUTINE lim_rhg *** … … 98 98 !! e.g. in the Canadian Archipelago 99 99 !! 100 !! ** References : Hunke and Dukowicz, JPO97 101 !! Bouillon et al., 08, in prep (update this when 102 !! published) 103 !! Vancoppenolle et al., OM08 100 !! References : Hunke and Dukowicz, JPO97 101 !! Bouillon et al., Ocean Modelling 2009 102 !! Vancoppenolle et al., Ocean Modelling 2008 103 !!------------------------------------------------------------------- 104 INTEGER, INTENT(in) :: k_j1 ! southern j-index for ice computation 105 INTEGER, INTENT(in) :: k_jpj ! northern j-index for ice computation 104 106 !! 105 !!------------------------------------------------------------------- 106 ! * Arguments 107 ! 108 INTEGER, INTENT(in) :: & 109 k_j1 , & !: southern j-index for ice computation 110 k_jpj !: northern j-index for ice computation 111 112 ! * Local variables 113 INTEGER :: ji, jj !: dummy loop indices 114 115 INTEGER :: & 116 jter !: temporary integers 117 107 INTEGER :: ji, jj ! dummy loop indices 108 INTEGER :: jter ! local integers 118 109 CHARACTER (len=50) :: charout 119 120 REAL(wp) :: & 121 zt11, zt12, zt21, zt22, & !: temporary scalars 122 ztagnx, ztagny, & !: wind stress on U/V points 123 delta ! 124 125 REAL(wp) :: & 126 za, & !: 127 zstms, & !: temporary scalar for ice strength 128 zsang, & !: temporary scalar for coriolis term 129 zmask !: mask for the computation of ice mass 110 REAL(wp) :: zt11, zt12, zt21, zt22, ztagnx, ztagny, delta ! 111 REAL(wp) :: za, zstms, zsang, zmask ! local scalars 130 112 131 113 REAL(wp),DIMENSION(jpi,jpj) :: & … … 145 127 146 128 REAL(wp) :: & 147 dtevp, & !: time step for subcycling 148 dtotel, & !: 149 ecc2, & !: square of yield ellipse eccenticity 150 z0, & !: temporary scalar 151 zr, & !: temporary scalar 152 zcca, zccb, & !: temporary scalars 153 zu_ice2, & !: 154 zv_ice1, & !: 155 zddc, zdtc, & !: temporary array for delta on corners 156 zdst, & !: temporary array for delta on centre 157 zdsshx, zdsshy, & !: term for the gradient of ocean surface 158 sigma1, sigma2 !: internal ice stress 129 dtevp, & ! time step for subcycling 130 dtotel, & ! 131 ecc2, & ! square of yield ellipse eccenticity 132 z0, & ! temporary scalar 133 zr, & ! temporary scalar 134 zcca, zccb, & ! temporary scalars 135 zu_ice2, & ! 136 zv_ice1, & ! 137 zddc, zdtc, & ! temporary array for delta on corners 138 zdst, & ! temporary array for delta on centre 139 zdsshx, zdsshy, & ! term for the gradient of ocean surface 140 sigma1, sigma2 ! internal ice stress 141 142 REAL(wp),DIMENSION(jpi,jpj) :: zf1, zf2 ! arrays for internal stresses 159 143 160 144 REAL(wp),DIMENSION(jpi,jpj) :: & 161 zf1, zf2 !: arrays for internal stresses 162 163 REAL(wp),DIMENSION(jpi,jpj) :: & 164 zdd, zdt, & !: Divergence and tension at centre of grid cells 165 zds, & !: Shear on northeast corner of grid cells 166 deltat, & !: Delta at centre of grid cells 167 deltac, & !: Delta on corners 168 zs1, zs2, & !: Diagonal stress tensor components zs1 and zs2 169 zs12 !: Non-diagonal stress tensor component zs12 145 zdd, zdt, & ! Divergence and tension at centre of grid cells 146 zds, & ! Shear on northeast corner of grid cells 147 deltat, & ! Delta at centre of grid cells 148 deltac, & ! Delta on corners 149 zs1, zs2, & ! Diagonal stress tensor components zs1 and zs2 150 zs12 ! Non-diagonal stress tensor component zs12 170 151 171 152 REAL(wp) :: & 172 zresm , & !: Maximal error on ice velocity 173 zindb , & !: ice (1) or not (0) 174 zdummy !: dummy argument 175 176 REAL(wp),DIMENSION(jpi,jpj) :: & 177 zu_ice , & !: Ice velocity on previous time step 178 zv_ice , & 179 zresr !: Local error on velocity 180 153 zresm , & ! Maximal error on ice velocity 154 zindb , & ! ice (1) or not (0) 155 zdummy ! dummy argument 156 157 REAL(wp),DIMENSION(jpi,jpj) :: zu_ice, zv_ice, zresr ! Local error on velocity 158 !!------------------------------------------------------------------- 159 #if defined key_lim2 && ! defined key_lim2_vp 160 # if defined key_agrif 161 USE ice_2, vt_s => hsnm 162 USE ice_2, vt_i => hicm 163 # else 164 vt_s => hsnm 165 vt_i => hicm 166 # endif 167 at_i(:,:) = 1. - frld(:,:) 168 #endif 181 169 ! 182 170 !------------------------------------------------------------------------------! … … 185 173 ! 186 174 ! Put every vector to 0 187 zpresh (:,:) = 0.0 ; zc1(:,:) = 0.0188 zpreshc(:,:) = 0. 0189 u_ice2 (:,:) = 0.0 ; v_ice1(:,:) = 0.0190 zdd (:,:) = 0.0 ; zdt(:,:) = 0.0 ; zds(:,:) = 0.0191 192 ! Ice strength on T-points 193 CALL lim_itd_me_icestrength( ridge_scheme_swi)194 195 ! Ice mass and temp variables 196 !CDIR NOVERRCHK 197 DO jj = k_j1 , k_jpj 175 zpresh (:,:) = 0._wp ; zc1 (:,:) = 0._wp 176 zpreshc(:,:) = 0._wp 177 u_ice2 (:,:) = 0._wp ; v_ice1(:,:) = 0._wp 178 zdd (:,:) = 0._wp ; zdt (:,:) = 0._wp ; zds(:,:) = 0._wp 179 180 #if defined key_lim3 181 CALL lim_itd_me_icestrength( ridge_scheme_swi ) ! LIM-3: Ice strength on T-points 182 #endif 183 184 !CDIR NOVERRCHK 185 DO jj = k_j1 , k_jpj ! Ice mass and temp variables 198 186 !CDIR NOVERRCHK 199 187 DO ji = 1 , jpi 200 188 zc1(ji,jj) = tms(ji,jj) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj) ) 201 zpresh(ji,jj) = tms(ji,jj) * strength(ji,jj) / 2. 189 #if defined key_lim3 190 zpresh(ji,jj) = tms(ji,jj) * strength(ji,jj) * 0.5_wp 191 #endif 202 192 ! tmi = 1 where there is ice or on land 203 tmi(ji,jj) = 1.0 - ( 1.0 - MAX( 0.0 , SIGN ( 1.0 , vt_i(ji,jj) - & 204 epsd ) ) ) * tms(ji,jj) 193 tmi(ji,jj) = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - epsd ) ) ) * tms(ji,jj) 205 194 END DO 206 195 END DO … … 247 236 DO ji = fs_2, fs_jpim1 248 237 249 zt11 = tms(ji ,jj)*e1t(ji,jj)250 zt12 = tms(ji+1,jj) *e1t(ji+1,jj)251 zt21 = tms(ji,jj )*e2t(ji,jj)252 zt22 = tms(ji,jj+1) *e2t(ji,jj+1)238 zt11 = tms(ji ,jj) * e1t(ji ,jj) 239 zt12 = tms(ji+1,jj) * e1t(ji+1,jj) 240 zt21 = tms(ji,jj ) * e2t(ji,jj ) 241 zt22 = tms(ji,jj+1) * e2t(ji,jj+1) 253 242 254 243 ! Leads area. 255 zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + & 256 & zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + epsd ) 257 zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + & 258 & zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + epsd ) 244 zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + epsd ) 245 zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + epsd ) 259 246 260 247 ! Mass, coriolis coeff. and currents 261 248 zmass1(ji,jj) = ( zt12*zc1(ji,jj) + zt11*zc1(ji+1,jj) ) / (zt11+zt12+epsd) 262 249 zmass2(ji,jj) = ( zt22*zc1(ji,jj) + zt21*zc1(ji,jj+1) ) / (zt21+zt22+epsd) 263 zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj)*fcor(ji,jj) + & 264 e1t(ji,jj)*fcor(ji+1,jj) ) & 265 / (e1t(ji,jj) + e1t(ji+1,jj) + epsd ) 266 zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1)*fcor(ji,jj) + & 267 e2t(ji,jj)*fcor(ji,jj+1) ) & 268 / ( e2t(ji,jj+1) + e2t(ji,jj) + epsd ) 250 zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj)*fcor(ji,jj) + e1t(ji,jj)*fcor(ji+1,jj) ) & 251 & / ( e1t(ji,jj) + e1t(ji+1,jj) + epsd ) 252 zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1)*fcor(ji,jj) + e2t(ji,jj)*fcor(ji,jj+1) ) & 253 & / ( e2t(ji,jj+1) + e2t(ji,jj) + epsd ) 269 254 ! 270 255 u_oce1(ji,jj) = u_oce(ji,jj) … … 290 275 ! include it later 291 276 292 zdsshx = ( ssh_m(ji+1,jj) - ssh_m(ji,jj))/e1u(ji,jj)293 zdsshy = ( ssh_m(ji,jj+1) - ssh_m(ji,jj))/e2v(ji,jj)277 zdsshx = ( ssh_m(ji+1,jj) - ssh_m(ji,jj) ) / e1u(ji,jj) 278 zdsshy = ( ssh_m(ji,jj+1) - ssh_m(ji,jj) ) / e2v(ji,jj) 294 279 295 280 za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx … … 306 291 ! Time step for subcycling 307 292 dtevp = rdt_ice / nevp 308 dtotel = dtevp / ( 2. 0* telast )293 dtotel = dtevp / ( 2._wp * telast ) 309 294 310 295 !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter) 311 ecc2 = ecc *ecc296 ecc2 = ecc * ecc 312 297 313 298 !-Initialise stress tensor 314 zs1 (:,:) = stress1_i(:,:)315 zs2 (:,:) = stress2_i(:,:)299 zs1 (:,:) = stress1_i (:,:) 300 zs2 (:,:) = stress2_i (:,:) 316 301 zs12(:,:) = stress12_i(:,:) 317 302 318 ! ----------------------!303 ! !----------------------! 319 304 DO jter = 1 , nevp ! loop over jter ! 320 ! ----------------------!305 ! !----------------------! 321 306 DO jj = k_j1, k_jpj-1 322 zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step307 zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 323 308 zv_ice(:,jj) = v_ice(:,jj) 324 309 END DO … … 385 370 END DO 386 371 END DO 387 CALL lbc_lnk( v_ice1(:,:), 'U', -1. ) 388 CALL lbc_lnk( u_ice2(:,:), 'V', -1. ) 372 CALL lbc_lnk( v_ice1, 'U', -1. ) ; CALL lbc_lnk( u_ice2, 'V', -1. ) ! lateral boundary cond. 389 373 390 374 !CDIR NOVERRCHK … … 394 378 395 379 !- Calculate Delta at centre of grid cells 396 zdst = ( e2u( ji , jj ) * v_ice1(ji,jj)&397 & - e2u( ji-1, jj) * v_ice1(ji-1,jj) &398 & + e1v( ji , jj ) * u_ice2(ji,jj)&399 & - e1v( ji , jj-1) * u_ice2(ji,jj-1) &400 & ) 380 zdst = ( e2u(ji , jj) * v_ice1(ji ,jj) & 381 & - e2u(ji-1, jj) * v_ice1(ji-1,jj) & 382 & + e1v(ji, jj ) * u_ice2(ji,jj ) & 383 & - e1v(ji, jj-1) * u_ice2(ji,jj-1) & 384 & ) & 401 385 & / area(ji,jj) 402 386 403 delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + & 404 & ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 ) 405 deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 + & 406 (zdt(ji,jj)**2 + zdst**2)*usecc2), creepl ) 387 delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 ) 388 deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 + (zdt(ji,jj)**2 + zdst**2)*usecc2), creepl ) 407 389 408 390 !-Calculate stress tensor components zs1 and zs2
Note: See TracChangeset
for help on using the changeset viewer.
