- Timestamp:
- 2015-07-08T17:13:59+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5107_hadgem3_mct/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
r5279 r5568 6 6 !! history : 1.0 ! 2002-08 (C. Ethe, G. Madec) original VP code 7 7 !! 3.0 ! 2007-03 (MA Morales Maqueda, S. Bouillon, M. Vancoppenolle) LIM3: EVP-Cgrid 8 !! 4.0! 2011-02 (G. Madec) dynamical allocation8 !! 3.5 ! 2011-02 (G. Madec) dynamical allocation 9 9 !!---------------------------------------------------------------------- 10 10 #if defined key_lim3 … … 20 20 USE sbc_ice ! Surface boundary condition: ice fields 21 21 USE ice ! LIM-3 variables 22 USE par_ice ! LIM-3 parameters23 22 USE dom_ice ! LIM-3 domain 24 23 USE limrhg ! LIM-3 rheology … … 31 30 USE timing ! Timing 32 31 USE limcons ! conservation tests 32 USE limvar 33 33 34 34 IMPLICIT NONE … … 76 76 CALL wrk_alloc( jpj, zswitch, zmsk ) 77 77 78 CALL lim_var_agg(1) ! aggregate ice categories 79 78 80 IF( kt == nit000 ) CALL lim_dyn_init ! Initialization (first time-step only) 79 81 … … 83 85 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 84 86 85 u_ice_b(:,:) = u_ice(:,:) * tmu(:,:)86 v_ice_b(:,:) = v_ice(:,:) * tmv(:,:)87 u_ice_b(:,:) = u_ice(:,:) * umask(:,:,1) 88 v_ice_b(:,:) = v_ice(:,:) * vmask(:,:,1) 87 89 88 90 ! Rheology (ice dynamics) … … 101 103 DO jj = 1, jpj 102 104 zswitch(jj) = SUM( 1.0 - at_i(:,jj) ) ! = REAL(jpj) if ocean everywhere on a j-line 103 zmsk (jj) = SUM( tmask(:,jj,1) ) ! = 0 if land everywhere on a j-line105 zmsk (jj) = SUM( tmask(:,jj,1) ) ! = 0 if land everywhere on a j-line 104 106 END DO 105 107 … … 157 159 zv_io(:,:) = v_ice(:,:) - ssv_m(:,:) 158 160 ! frictional velocity at T-point 159 zcoef = 0.5_wp * cw161 zcoef = 0.5_wp * rn_cio 160 162 DO jj = 2, jpjm1 161 163 DO ji = fs_2, fs_jpim1 ! vector opt. 162 164 ust2s(ji,jj) = zcoef * ( zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj) & 163 & + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) * tms(ji,jj)165 & + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) * tmask(ji,jj,1) 164 166 END DO 165 167 END DO … … 170 172 ELSE ! no ice dynamics : transmit directly the atmospheric stress to the ocean 171 173 ! 172 zcoef = SQRT( 0.5_wp ) /rau0174 zcoef = SQRT( 0.5_wp ) * r1_rau0 173 175 DO jj = 2, jpjm1 174 176 DO ji = fs_2, fs_jpim1 ! vector opt. 175 177 ust2s(ji,jj) = zcoef * SQRT( utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj) & 176 & + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) * tms(ji,jj)178 & + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) * tmask(ji,jj,1) 177 179 END DO 178 180 END DO … … 189 191 CALL prt_ctl(tab2d_1=delta_i , clinfo1=' lim_dyn : delta_i :') 190 192 CALL prt_ctl(tab2d_1=strength , clinfo1=' lim_dyn : strength :') 191 CALL prt_ctl(tab2d_1= area, clinfo1=' lim_dyn : cell area :')193 CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_dyn : cell area :') 192 194 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_dyn : at_i :') 193 195 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_dyn : vt_i :') … … 241 243 !!------------------------------------------------------------------- 242 244 INTEGER :: ios ! Local integer output status for namelist read 243 NAMELIST/namicedyn/ epsd, om, cw, pstar, & 244 & c_rhg, creepl, ecc, ahi0, & 245 & nevp, relast, alphaevp, hminrhg 245 NAMELIST/namicedyn/ nn_icestr, ln_icestr_bvf, rn_pe_rdg, rn_pstar, rn_crhg, rn_cio, rn_creepl, rn_ecc, & 246 & nn_nevp, rn_relast, nn_ahi0, rn_ahi0_ref 247 INTEGER :: ji, jj 248 REAL(wp) :: za00, zd_max 246 249 !!------------------------------------------------------------------- 247 250 … … 259 262 WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics ' 260 263 WRITE(numout,*) '~~~~~~~~~~~~' 261 WRITE(numout,*) ' tolerance parameter epsd = ', epsd262 WRITE(numout,*) ' relaxation constant om = ', om263 WRITE(numout,*) ' drag coefficient for oceanic stress cw = ', cw264 WRITE(numout,*) ' first bulk-rheology parameter pstar = ', pstar265 WRITE(numout,*) ' second bulk-rhelogy parameter c_rhg = ', c_rhg266 WRITE(numout,*) ' creep limit creepl = ', creepl267 WRITE(numout,*) ' eccentricity of the elliptical yield curve ecc = ', ecc268 WRITE(numout,*) ' horizontal diffusivity coeff. for sea-ice ahi0 = ', ahi0269 WRITE(numout,*) ' number of iterations for subcycling nevp = ',nevp270 WRITE(numout,*) ' ratio of elastic timescale over ice time step relast = ',relast271 WRITE(numout,*) ' coefficient for the solution of int. stresses alphaevp = ', alphaevp272 WRITE(numout,*) ' min ice thickness for rheology calculations hminrhg = ', hminrhg264 WRITE(numout,*)' ice strength parameterization (0=Hibler 1=Rothrock) nn_icestr = ', nn_icestr 265 WRITE(numout,*)' Including brine volume in ice strength comp. ln_icestr_bvf = ', ln_icestr_bvf 266 WRITE(numout,*)' Ratio of ridging work to PotEner change in ridging rn_pe_rdg = ', rn_pe_rdg 267 WRITE(numout,*) ' drag coefficient for oceanic stress rn_cio = ', rn_cio 268 WRITE(numout,*) ' first bulk-rheology parameter rn_pstar = ', rn_pstar 269 WRITE(numout,*) ' second bulk-rhelogy parameter rn_crhg = ', rn_crhg 270 WRITE(numout,*) ' creep limit rn_creepl = ', rn_creepl 271 WRITE(numout,*) ' eccentricity of the elliptical yield curve rn_ecc = ', rn_ecc 272 WRITE(numout,*) ' number of iterations for subcycling nn_nevp = ', nn_nevp 273 WRITE(numout,*) ' ratio of elastic timescale over ice time step rn_relast = ', rn_relast 274 WRITE(numout,*) ' horizontal diffusivity calculation nn_ahi0 = ', nn_ahi0 275 WRITE(numout,*) ' horizontal diffusivity coeff. (orca2 grid) rn_ahi0_ref = ', rn_ahi0_ref 273 276 ENDIF 274 277 ! 275 usecc2 = 1._wp / ( ecc * ecc ) 276 rhoco = rau0 * cw 277 278 ! elastic damping 279 telast = relast * rdt_ice 280 281 ! Diffusion coefficients. 282 ahiu(:,:) = ahi0 * umask(:,:,1) 283 ahiv(:,:) = ahi0 * vmask(:,:,1) 284 ! 278 usecc2 = 1._wp / ( rn_ecc * rn_ecc ) 279 rhoco = rau0 * rn_cio 280 ! 281 ! Diffusion coefficients 282 SELECT CASE( nn_ahi0 ) 283 284 CASE( 0 ) 285 ahiu(:,:) = rn_ahi0_ref 286 ahiv(:,:) = rn_ahi0_ref 287 288 IF(lwp) WRITE(numout,*) '' 289 IF(lwp) WRITE(numout,*) ' laplacian operator: ahim constant = rn_ahi0_ref' 290 291 CASE( 1 ) 292 293 zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 294 IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain 295 296 ahiu(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp ! 1.e05 = 100km = max grid space at 60° latitude in orca2 297 ! (60° = min latitude for ice cover) 298 ahiv(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp 299 300 IF(lwp) WRITE(numout,*) '' 301 IF(lwp) WRITE(numout,*) ' laplacian operator: ahim proportional to max of e1 e2 over the domain (', zd_max, ')' 302 IF(lwp) WRITE(numout,*) ' value for ahim = ', rn_ahi0_ref * zd_max * 1.e-05_wp 303 304 CASE( 2 ) 305 306 zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 307 IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain 308 309 za00 = rn_ahi0_ref * 1.e-05_wp ! 1.e05 = 100km = max grid space at 60° latitude in orca2 310 ! (60° = min latitude for ice cover) 311 DO jj = 1, jpj 312 DO ji = 1, jpi 313 ahiu(ji,jj) = za00 * MAX( e1t(ji,jj), e2t(ji,jj) ) * umask(ji,jj,1) 314 ahiv(ji,jj) = za00 * MAX( e1f(ji,jj), e2f(ji,jj) ) * vmask(ji,jj,1) 315 END DO 316 END DO 317 ! 318 IF(lwp) WRITE(numout,*) '' 319 IF(lwp) WRITE(numout,*) ' laplacian operator: ahim proportional to e1' 320 IF(lwp) WRITE(numout,*) ' maximum grid-spacing = ', zd_max, ' maximum value for ahim = ', za00*zd_max 321 322 END SELECT 323 285 324 END SUBROUTINE lim_dyn_init 286 325
Note: See TracChangeset
for help on using the changeset viewer.
