- Timestamp:
- 2015-01-20T15:26:13+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
r4161 r5038 30 30 USE lib_fortran ! glob_sum 31 31 USE timing ! Timing 32 USE limcons ! conservation tests 32 33 33 34 IMPLICIT NONE … … 63 64 INTEGER :: i_j1, i_jpj ! Starting/ending j-indices for rheology 64 65 REAL(wp) :: zcoef ! local scalar 65 REAL(wp), POINTER, DIMENSION(:) :: z ind! i-averaged indicator of sea-ice66 REAL(wp), POINTER, DIMENSION(:) :: zswitch ! i-averaged indicator of sea-ice 66 67 REAL(wp), POINTER, DIMENSION(:) :: zmsk ! i-averaged of tmask 67 68 REAL(wp), POINTER, DIMENSION(:,:) :: zu_io, zv_io ! ice-ocean velocity 68 REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset)69 REAL(wp) :: z chk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset)69 ! 70 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 70 71 !!--------------------------------------------------------------------- 71 72 … … 73 74 74 75 CALL wrk_alloc( jpi, jpj, zu_io, zv_io ) 75 CALL wrk_alloc( jpj, zind, zmsk ) 76 77 ! ------------------------------- 78 !- check conservation (C Rousset) 79 IF (ln_limdiahsb) THEN 80 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 81 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 82 zchk_fw_b = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 83 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 84 ENDIF 85 !- check conservation (C Rousset) 86 ! ------------------------------- 76 CALL wrk_alloc( jpj, zswitch, zmsk ) 87 77 88 78 IF( kt == nit000 ) CALL lim_dyn_init ! Initialization (first time-step only) … … 90 80 IF( ln_limdyn ) THEN 91 81 ! 92 old_u_ice(:,:) = u_ice(:,:) * tmu(:,:) 93 old_v_ice(:,:) = v_ice(:,:) * tmv(:,:) 82 ! conservation test 83 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 84 85 u_ice_b(:,:) = u_ice(:,:) * tmu(:,:) 86 v_ice_b(:,:) = v_ice(:,:) * tmv(:,:) 94 87 95 88 ! Rheology (ice dynamics) … … 107 100 ! 108 101 DO jj = 1, jpj 109 z ind(jj) = SUM( 1.0 - at_i(:,jj) ) ! = REAL(jpj) if ocean everywhere on a j-line102 zswitch(jj) = SUM( 1.0 - at_i(:,jj) ) ! = REAL(jpj) if ocean everywhere on a j-line 110 103 zmsk(jj) = SUM( tmask(:,jj,1) ) ! = 0 if land everywhere on a j-line 111 104 END DO … … 117 110 i_j1 = njeq 118 111 i_jpj = jpj 119 DO WHILE ( i_j1 <= jpj .AND. z ind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )112 DO WHILE ( i_j1 <= jpj .AND. zswitch(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 ) 120 113 i_j1 = i_j1 + 1 121 114 END DO … … 127 120 i_j1 = 1 128 121 i_jpj = njeq 129 DO WHILE ( i_jpj >= 1 .AND. z ind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )122 DO WHILE ( i_jpj >= 1 .AND. zswitch(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 ) 130 123 i_jpj = i_jpj - 1 131 124 END DO … … 139 132 ! ! latitude strip 140 133 i_j1 = 1 141 DO WHILE ( i_j1 <= jpj .AND. z ind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )134 DO WHILE ( i_j1 <= jpj .AND. zswitch(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 ) 142 135 i_j1 = i_j1 + 1 143 136 END DO … … 145 138 146 139 i_jpj = jpj 147 DO WHILE ( i_jpj >= 1 .AND. z ind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )140 DO WHILE ( i_jpj >= 1 .AND. zswitch(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 ) 148 141 i_jpj = i_jpj - 1 149 142 END DO … … 171 164 END DO 172 165 END DO 166 ! 167 ! conservation test 168 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 173 169 ! 174 170 ELSE ! no ice dynamics : transmit directly the atmospheric stress to the ocean … … 224 220 ENDIF 225 221 ! 226 ! -------------------------------227 !- check conservation (C Rousset)228 IF (ln_limdiahsb) THEN229 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b230 zchk_fw = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b231 232 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) / rdt_ice233 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic )234 235 zchk_vmin = glob_min(v_i)236 zchk_amax = glob_max(SUM(a_i,dim=3))237 zchk_amin = glob_min(a_i)238 239 IF(lwp) THEN240 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limdyn) = ',(zchk_v_i * rday)241 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limdyn) = ',(zchk_smv * rday)242 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limdyn) = ',(zchk_vmin * 1.e-3)243 !IF ( zchk_amax > amax+1.e-10 ) WRITE(numout,*) 'violation a_i>amax (limdyn) = ',zchk_amax244 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limdyn) = ',zchk_amin245 ENDIF246 ENDIF247 !- check conservation (C Rousset)248 ! -------------------------------249 250 222 CALL wrk_dealloc( jpi, jpj, zu_io, zv_io ) 251 CALL wrk_dealloc( jpj, z ind, zmsk )223 CALL wrk_dealloc( jpj, zswitch, zmsk ) 252 224 ! 253 225 IF( nn_timing == 1 ) CALL timing_stop('limdyn') … … 269 241 !!------------------------------------------------------------------- 270 242 INTEGER :: ios ! Local integer output status for namelist read 271 NAMELIST/namicedyn/ epsd, alpha, & 272 & dm, nbiter, nbitdr, om, resl, cw, angvg, pstar, & 273 & c_rhg, etamn, creepl, ecc, ahi0, & 274 & nevp, telast, alphaevp, hminrhg 243 NAMELIST/namicedyn/ epsd, om, cw, pstar, & 244 & c_rhg, creepl, ecc, ahi0, & 245 & nevp, relast, alphaevp, hminrhg 275 246 !!------------------------------------------------------------------- 276 247 … … 282 253 READ ( numnam_ice_cfg, namicedyn, IOSTAT = ios, ERR = 902 ) 283 254 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicedyn in configuration namelist', lwp ) 284 WRITE ( numoni, namicedyn )255 IF(lwm) WRITE ( numoni, namicedyn ) 285 256 286 257 IF(lwp) THEN ! control print … … 289 260 WRITE(numout,*) '~~~~~~~~~~~~' 290 261 WRITE(numout,*) ' tolerance parameter epsd = ', epsd 291 WRITE(numout,*) ' coefficient for semi-implicit coriolis alpha = ', alpha292 WRITE(numout,*) ' diffusion constant for dynamics dm = ', dm293 WRITE(numout,*) ' number of sub-time steps for relaxation nbiter = ', nbiter294 WRITE(numout,*) ' maximum number of iterations for relaxation nbitdr = ', nbitdr295 262 WRITE(numout,*) ' relaxation constant om = ', om 296 WRITE(numout,*) ' maximum value for the residual of relaxation resl = ', resl297 263 WRITE(numout,*) ' drag coefficient for oceanic stress cw = ', cw 298 WRITE(numout,*) ' turning angle for oceanic stress angvg = ', angvg299 264 WRITE(numout,*) ' first bulk-rheology parameter pstar = ', pstar 300 265 WRITE(numout,*) ' second bulk-rhelogy parameter c_rhg = ', c_rhg 301 WRITE(numout,*) ' minimun value for viscosity etamn = ', etamn302 266 WRITE(numout,*) ' creep limit creepl = ', creepl 303 267 WRITE(numout,*) ' eccentricity of the elliptical yield curve ecc = ', ecc 304 268 WRITE(numout,*) ' horizontal diffusivity coeff. for sea-ice ahi0 = ', ahi0 305 269 WRITE(numout,*) ' number of iterations for subcycling nevp = ', nevp 306 WRITE(numout,*) ' timescale for elastic waves telast = ', telast270 WRITE(numout,*) ' ratio of elastic timescale over ice time step relast = ', relast 307 271 WRITE(numout,*) ' coefficient for the solution of int. stresses alphaevp = ', alphaevp 308 272 WRITE(numout,*) ' min ice thickness for rheology calculations hminrhg = ', hminrhg 309 273 ENDIF 310 274 ! 311 IF( angvg /= 0._wp ) THEN312 CALL ctl_warn( 'lim_dyn_init: turning angle for oceanic stress not properly coded for EVP ', &313 & '(see limsbc module). We force angvg = 0._wp' )314 angvg = 0._wp315 ENDIF316 317 275 usecc2 = 1._wp / ( ecc * ecc ) 318 276 rhoco = rau0 * cw 319 angvg = angvg * rad 320 sangvg = SIN( angvg ) 321 cangvg = COS( angvg ) 322 pstarh = pstar * 0.5_wp 277 278 ! elastic damping 279 telast = relast * rdt_ice 323 280 324 281 ! Diffusion coefficients.
Note: See TracChangeset
for help on using the changeset viewer.
