- 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/limdyn.F90
r4688 r6225 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 … … 64 64 INTEGER :: i_j1, i_jpj ! Starting/ending j-indices for rheology 65 65 REAL(wp) :: zcoef ! local scalar 66 REAL(wp), POINTER, DIMENSION(:) :: z ind! i-averaged indicator of sea-ice66 REAL(wp), POINTER, DIMENSION(:) :: zswitch ! i-averaged indicator of sea-ice 67 67 REAL(wp), POINTER, DIMENSION(:) :: zmsk ! i-averaged of tmask 68 68 REAL(wp), POINTER, DIMENSION(:,:) :: zu_io, zv_io ! ice-ocean velocity … … 74 74 75 75 CALL wrk_alloc( jpi, jpj, zu_io, zv_io ) 76 CALL wrk_alloc( jpj, zind, zmsk ) 76 CALL wrk_alloc( jpj, zswitch, zmsk ) 77 78 CALL lim_var_agg(1) ! aggregate ice categories 77 79 78 80 IF( kt == nit000 ) CALL lim_dyn_init ! Initialization (first time-step only) … … 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 old_u_ice(:,:) = u_ice(:,:) * tmu(:,:)86 old_v_ice(:,:) = 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) … … 100 102 ! 101 103 DO jj = 1, jpj 102 z ind(jj) = SUM( 1.0 - at_i(:,jj) ) ! = REAL(jpj) if ocean everywhere on a j-line103 zmsk (jj) = SUM( tmask(:,jj,1) ) ! = 0 if land everywhere on a j-line104 zswitch(jj) = SUM( 1.0 - at_i(:,jj) ) ! = REAL(jpj) if ocean everywhere on a j-line 105 zmsk (jj) = SUM( tmask(:,jj,1) ) ! = 0 if land everywhere on a j-line 104 106 END DO 105 107 … … 110 112 i_j1 = njeq 111 113 i_jpj = jpj 112 DO WHILE ( i_j1 <= jpj .AND. z ind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )114 DO WHILE ( i_j1 <= jpj .AND. zswitch(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 ) 113 115 i_j1 = i_j1 + 1 114 116 END DO … … 120 122 i_j1 = 1 121 123 i_jpj = njeq 122 DO WHILE ( i_jpj >= 1 .AND. z ind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )124 DO WHILE ( i_jpj >= 1 .AND. zswitch(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 ) 123 125 i_jpj = i_jpj - 1 124 126 END DO … … 132 134 ! ! latitude strip 133 135 i_j1 = 1 134 DO WHILE ( i_j1 <= jpj .AND. z ind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )136 DO WHILE ( i_j1 <= jpj .AND. zswitch(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 ) 135 137 i_j1 = i_j1 + 1 136 138 END DO … … 138 140 139 141 i_jpj = jpj 140 DO WHILE ( i_jpj >= 1 .AND. z ind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )142 DO WHILE ( i_jpj >= 1 .AND. zswitch(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 ) 141 143 i_jpj = i_jpj - 1 142 144 END DO … … 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=e1e2t , 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 :') … … 221 223 ! 222 224 CALL wrk_dealloc( jpi, jpj, zu_io, zv_io ) 223 CALL wrk_dealloc( jpj, z ind, zmsk )225 CALL wrk_dealloc( jpj, zswitch, zmsk ) 224 226 ! 225 227 IF( nn_timing == 1 ) CALL timing_stop('limdyn') … … 241 243 !!------------------------------------------------------------------- 242 244 INTEGER :: ios ! Local integer output status for namelist read 243 NAMELIST/namicedyn/ epsd, om, cw, angvg, pstar, & 244 & c_rhg, creepl, ecc, ahi0, & 245 & nevp, telast, 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 = ', epsd 262 WRITE(numout,*) ' relaxation constant om = ', om 263 WRITE(numout,*) ' drag coefficient for oceanic stress cw = ', cw 264 WRITE(numout,*) ' turning angle for oceanic stress angvg = ', angvg 265 WRITE(numout,*) ' first bulk-rheology parameter pstar = ', pstar 266 WRITE(numout,*) ' second bulk-rhelogy parameter c_rhg = ', c_rhg 267 WRITE(numout,*) ' creep limit creepl = ', creepl 268 WRITE(numout,*) ' eccentricity of the elliptical yield curve ecc = ', ecc 269 WRITE(numout,*) ' horizontal diffusivity coeff. for sea-ice ahi0 = ', ahi0 270 WRITE(numout,*) ' number of iterations for subcycling nevp = ', nevp 271 WRITE(numout,*) ' timescale for elastic waves telast = ', telast 272 WRITE(numout,*) ' coefficient for the solution of int. stresses alphaevp = ', alphaevp 273 WRITE(numout,*) ' min ice thickness for rheology calculations hminrhg = ', hminrhg 264 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 274 276 ENDIF 275 277 ! 276 IF( angvg /= 0._wp ) THEN 277 CALL ctl_warn( 'lim_dyn_init: turning angle for oceanic stress not properly coded for EVP ', & 278 & '(see limsbc module). We force angvg = 0._wp' ) 279 angvg = 0._wp 280 ENDIF 281 282 usecc2 = 1._wp / ( ecc * ecc ) 283 rhoco = rau0 * cw 284 angvg = angvg * rad 285 sangvg = SIN( angvg ) 286 cangvg = COS( angvg ) 287 pstarh = pstar * 0.5_wp 288 289 ! Diffusion coefficients. 290 ahiu(:,:) = ahi0 * umask(:,:,1) 291 ahiv(:,:) = ahi0 * vmask(:,:,1) 292 ! 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 293 324 END SUBROUTINE lim_dyn_init 294 325
Note: See TracChangeset
for help on using the changeset viewer.