Changeset 5260 for branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
- Timestamp:
- 2015-05-12T12:37:15+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
r4624 r5260 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 … … 30 29 USE lib_fortran ! glob_sum 31 30 USE timing ! Timing 31 USE limcons ! conservation tests 32 USE limvar 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 ) 77 78 CALL lim_var_agg(1) ! aggregate ice categories 87 79 88 80 IF( kt == nit000 ) CALL lim_dyn_init ! Initialization (first time-step only) … … 90 82 IF( ln_limdyn ) THEN 91 83 ! 92 old_u_ice(:,:) = u_ice(:,:) * tmu(:,:) 93 old_v_ice(:,:) = v_ice(:,:) * tmv(:,:) 84 ! conservation test 85 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 86 87 u_ice_b(:,:) = u_ice(:,:) * umask(:,:,1) 88 v_ice_b(:,:) = v_ice(:,:) * vmask(:,:,1) 94 89 95 90 ! Rheology (ice dynamics) … … 107 102 ! 108 103 DO jj = 1, jpj 109 z ind(jj) = SUM( 1.0 - at_i(:,jj) ) ! = REAL(jpj) if ocean everywhere on a j-line110 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 111 106 END DO 112 107 … … 117 112 i_j1 = njeq 118 113 i_jpj = jpj 119 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 ) 120 115 i_j1 = i_j1 + 1 121 116 END DO … … 127 122 i_j1 = 1 128 123 i_jpj = njeq 129 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 ) 130 125 i_jpj = i_jpj - 1 131 126 END DO … … 139 134 ! ! latitude strip 140 135 i_j1 = 1 141 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 ) 142 137 i_j1 = i_j1 + 1 143 138 END DO … … 145 140 146 141 i_jpj = jpj 147 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 ) 148 143 i_jpj = i_jpj - 1 149 144 END DO … … 164 159 zv_io(:,:) = v_ice(:,:) - ssv_m(:,:) 165 160 ! frictional velocity at T-point 166 zcoef = 0.5_wp * cw161 zcoef = 0.5_wp * rn_cio 167 162 DO jj = 2, jpjm1 168 163 DO ji = fs_2, fs_jpim1 ! vector opt. 169 164 ust2s(ji,jj) = zcoef * ( zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj) & 170 & + 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) 171 166 END DO 172 167 END DO 173 168 ! 169 ! conservation test 170 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 171 ! 174 172 ELSE ! no ice dynamics : transmit directly the atmospheric stress to the ocean 175 173 ! 176 zcoef = SQRT( 0.5_wp ) /rau0174 zcoef = SQRT( 0.5_wp ) * r1_rau0 177 175 DO jj = 2, jpjm1 178 176 DO ji = fs_2, fs_jpim1 ! vector opt. 179 177 ust2s(ji,jj) = zcoef * SQRT( utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj) & 180 & + 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) 181 179 END DO 182 180 END DO … … 193 191 CALL prt_ctl(tab2d_1=delta_i , clinfo1=' lim_dyn : delta_i :') 194 192 CALL prt_ctl(tab2d_1=strength , clinfo1=' lim_dyn : strength :') 195 CALL prt_ctl(tab2d_1= area, clinfo1=' lim_dyn : cell area :')193 CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_dyn : cell area :') 196 194 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_dyn : at_i :') 197 195 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_dyn : vt_i :') … … 224 222 ENDIF 225 223 ! 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 224 CALL wrk_dealloc( jpi, jpj, zu_io, zv_io ) 251 CALL wrk_dealloc( jpj, z ind, zmsk )225 CALL wrk_dealloc( jpj, zswitch, zmsk ) 252 226 ! 253 227 IF( nn_timing == 1 ) CALL timing_stop('limdyn') … … 269 243 !!------------------------------------------------------------------- 270 244 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, hminrhg245 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 275 249 !!------------------------------------------------------------------- 276 250 … … 288 262 WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics ' 289 263 WRITE(numout,*) '~~~~~~~~~~~~' 290 WRITE(numout,*) ' tolerance parameter epsd = ', epsd 291 WRITE(numout,*) ' coefficient for semi-implicit coriolis alpha = ', alpha 292 WRITE(numout,*) ' diffusion constant for dynamics dm = ', dm 293 WRITE(numout,*) ' number of sub-time steps for relaxation nbiter = ', nbiter 294 WRITE(numout,*) ' maximum number of iterations for relaxation nbitdr = ', nbitdr 295 WRITE(numout,*) ' relaxation constant om = ', om 296 WRITE(numout,*) ' maximum value for the residual of relaxation resl = ', resl 297 WRITE(numout,*) ' drag coefficient for oceanic stress cw = ', cw 298 WRITE(numout,*) ' turning angle for oceanic stress angvg = ', angvg 299 WRITE(numout,*) ' first bulk-rheology parameter pstar = ', pstar 300 WRITE(numout,*) ' second bulk-rhelogy parameter c_rhg = ', c_rhg 301 WRITE(numout,*) ' minimun value for viscosity etamn = ', etamn 302 WRITE(numout,*) ' creep limit creepl = ', creepl 303 WRITE(numout,*) ' eccentricity of the elliptical yield curve ecc = ', ecc 304 WRITE(numout,*) ' horizontal diffusivity coeff. for sea-ice ahi0 = ', ahi0 305 WRITE(numout,*) ' number of iterations for subcycling nevp = ', nevp 306 WRITE(numout,*) ' timescale for elastic waves telast = ', telast 307 WRITE(numout,*) ' coefficient for the solution of int. stresses alphaevp = ', alphaevp 308 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 309 276 ENDIF 310 277 ! 311 IF( angvg /= 0._wp ) THEN 312 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._wp 315 ENDIF 316 317 usecc2 = 1._wp / ( ecc * ecc ) 318 rhoco = rau0 * cw 319 angvg = angvg * rad 320 sangvg = SIN( angvg ) 321 cangvg = COS( angvg ) 322 pstarh = pstar * 0.5_wp 323 324 ! Diffusion coefficients. 325 ahiu(:,:) = ahi0 * umask(:,:,1) 326 ahiv(:,:) = ahi0 * vmask(:,:,1) 327 ! 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 328 324 END SUBROUTINE lim_dyn_init 329 325
Note: See TracChangeset
for help on using the changeset viewer.
