- Timestamp:
- 2018-12-17T12:25:09+01:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn.F90
r10312 r10399 21 21 USE icecor ! sea-ice: corrections 22 22 USE icevar ! sea-ice: operations 23 USE icectl ! sea-ice: control prints 23 24 ! 24 25 USE in_out_manager ! I/O manager … … 37 38 INTEGER :: nice_dyn ! choice of the type of dynamics 38 39 ! ! associated indices: 39 INTEGER, PARAMETER :: np_dyn FULL= 1 ! full ice dynamics (rheology + advection + ridging/rafting + correction)40 INTEGER, PARAMETER :: np_dynALL = 1 ! full ice dynamics (rheology + advection + ridging/rafting + correction) 40 41 INTEGER, PARAMETER :: np_dynRHGADV = 2 ! pure dynamics (rheology + advection) 41 INTEGER, PARAMETER :: np_dynADV = 3 ! only advection w prescribed vel.(rn_uvice + advection) 42 INTEGER, PARAMETER :: np_dynADV1D = 3 ! only advection 1D - test case from Schar & Smolarkiewicz 1996 43 INTEGER, PARAMETER :: np_dynADV2D = 4 ! only advection 2D w prescribed vel.(rn_uvice + advection) 42 44 ! 43 45 ! ** namelist (namdyn) ** 44 LOGICAL :: ln_dynFULL ! full ice dynamics (rheology + advection + ridging/rafting + correction) 45 LOGICAL :: ln_dynRHGADV ! no ridge/raft & no corrections (rheology + advection) 46 LOGICAL :: ln_dynADV ! only advection w prescribed vel.(rn_uvice + advection) 47 REAL(wp) :: rn_uice ! prescribed u-vel (case np_dynADV) 48 REAL(wp) :: rn_vice ! prescribed v-vel (case np_dynADV) 46 LOGICAL :: ln_dynALL ! full ice dynamics (rheology + advection + ridging/rafting + correction) 47 LOGICAL :: ln_dynRHGADV ! no ridge/raft & no corrections (rheology + advection) 48 LOGICAL :: ln_dynADV1D ! only advection in 1D w ice convergence (test case from Schar & Smolarkiewicz 1996) 49 LOGICAL :: ln_dynADV2D ! only advection in 2D w prescribed vel. (rn_uvice + advection) 50 REAL(wp) :: rn_uice ! prescribed u-vel (case np_dynADV1D & np_dynADV2D) 51 REAL(wp) :: rn_vice ! prescribed v-vel (case np_dynADV1D & np_dynADV2D) 49 52 50 53 !! * Substitutions … … 73 76 INTEGER :: ji, jj, jl ! dummy loop indices 74 77 REAL(wp) :: zcoefu, zcoefv 75 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zh max78 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhi_max, zhs_max 76 79 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdivu_i 77 80 !!-------------------------------------------------------------------- 78 81 ! 79 IF( ln_timing ) CALL timing_start('icedyn') 82 ! controls 83 IF( ln_timing ) CALL timing_start('icedyn') ! timing 84 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icedyn', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 80 85 ! 81 86 IF( kt == nit000 .AND. lwp ) THEN … … 88 93 tau_icebfr(:,:) = 0._wp 89 94 DO jl = 1, jpl 90 WHERE( h_i (:,:,jl) > ht_n(:,:) * rn_depfra ) tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr91 END DO 92 ENDIF 93 ! 94 zhmax(:,:,:) = h_i_b(:,:,:)!-- Record max of the surrounding 9-pts ice thick. (for CALL Hbig)95 WHERE( h_i_b(:,:,jl) > ht_n(:,:) * rn_depfra ) tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 96 END DO 97 ENDIF 98 ! 99 ! !-- Record max of the surrounding 9-pts ice thick. (for CALL Hbig) 95 100 DO jl = 1, jpl 96 101 DO jj = 2, jpjm1 97 DO ji = 2, jpim1 98 !!gm use of MAXVAL here is very probably less efficient than expending the 9 values 99 zhmax(ji,jj,jl) = MAX( epsi20, MAXVAL( h_i_b(ji-1:ji+1,jj-1:jj+1,jl) ) ) 102 DO ji = fs_2, fs_jpim1 103 zhi_max(ji,jj,jl) = MAX( epsi20, h_i_b(ji,jj,jl), h_i_b(ji+1,jj ,jl), h_i_b(ji ,jj+1,jl), & 104 & h_i_b(ji-1,jj ,jl), h_i_b(ji ,jj-1,jl), & 105 & h_i_b(ji+1,jj+1,jl), h_i_b(ji-1,jj-1,jl), & 106 & h_i_b(ji+1,jj-1,jl), h_i_b(ji-1,jj+1,jl) ) 107 zhs_max(ji,jj,jl) = MAX( epsi20, h_s_b(ji,jj,jl), h_s_b(ji+1,jj ,jl), h_s_b(ji ,jj+1,jl), & 108 & h_s_b(ji-1,jj ,jl), h_s_b(ji ,jj-1,jl), & 109 & h_s_b(ji+1,jj+1,jl), h_s_b(ji-1,jj-1,jl), & 110 & h_s_b(ji+1,jj-1,jl), h_s_b(ji-1,jj+1,jl) ) 100 111 END DO 101 112 END DO 102 113 END DO 103 CALL lbc_lnk ( zhmax(:,:,:), 'T', 1. )114 CALL lbc_lnk_multi( zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1. ) 104 115 ! 105 116 ! 106 117 SELECT CASE( nice_dyn ) !-- Set which dynamics is running 107 118 108 CASE ( np_dyn FULL )!== all dynamical processes ==!109 CALL ice_dyn_rhg ( kt ) ! -- rheology110 CALL ice_dyn_adv ( kt ) ; CALL Hbig( zh max ) ! -- advection of ice + correction on ice thickness111 CALL ice_dyn_rdgrft( kt ) ! -- ridging/rafting112 CALL ice_cor ( kt , 1 ) ! -- Corrections119 CASE ( np_dynALL ) !== all dynamical processes ==! 120 CALL ice_dyn_rhg ( kt ) ! -- rheology 121 CALL ice_dyn_adv ( kt ) ; CALL Hbig( zhi_max, zhs_max ) ! -- advection of ice + correction on ice thickness 122 CALL ice_dyn_rdgrft( kt ) ! -- ridging/rafting 123 CALL ice_cor ( kt , 1 ) ! -- Corrections 113 124 114 125 CASE ( np_dynRHGADV ) !== no ridge/raft & no corrections ==! 115 CALL ice_dyn_rhg ( kt ) ! -- rheology116 CALL ice_dyn_adv ( kt ) ! -- advection of ice117 CALL Hpiling ! -- simple pile-up (replaces ridging/rafting)118 119 CASE ( np_dynADV ) !== pure advection ==! (prescribed velocities)126 CALL ice_dyn_rhg ( kt ) ! -- rheology 127 CALL ice_dyn_adv ( kt ) ; CALL Hbig( zhi_max, zhs_max ) ! -- advection of ice + correction on ice thickness 128 CALL Hpiling ! -- simple pile-up (replaces ridging/rafting) 129 130 CASE ( np_dynADV1D ) !== pure advection ==! (1D) 120 131 ALLOCATE( zdivu_i(jpi,jpj) ) 121 122 u_ice(:,:) = rn_uice * umask(:,:,1) 123 v_ice(:,:) = rn_vice * vmask(:,:,1) 124 !CALL RANDOM_NUMBER(u_ice(:,:)) ; u_ice(:,:) = u_ice(:,:) * 0.1 + rn_uice * 0.9 * umask(:,:,1) 125 !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1) 126 ! --- monotonicity test from Lipscomb et al 2004 --- ! 132 ! --- monotonicity test from Schar & Smolarkiewicz 1996 --- ! 127 133 ! CFL = 0.5 at a distance from the bound of 1/6 of the basin length 128 134 ! Then for dx = 2m and dt = 1s => rn_uice = u (1/6th) = 1m/s 129 !DO jj = 1, jpj130 !DO ji = 1, jpi131 !zcoefu = ( REAL(jpiglo+1)*0.5 - REAL(ji+nimpp-1) ) / ( REAL(jpiglo+1)*0.5 - 1. )132 !zcoefv = ( REAL(jpjglo+1)*0.5 - REAL(jj+njmpp-1) ) / ( REAL(jpjglo+1)*0.5 - 1. )133 !u_ice(ji,jj) = rn_uice * 1.5 * SIGN( 1., zcoefu ) * ABS( zcoefu ) * umask(ji,jj,1)134 !v_ice(ji,jj) = rn_vice * 1.5 * SIGN( 1., zcoefv ) * ABS( zcoefv ) * vmask(ji,jj,1)135 !END DO136 !END DO135 DO jj = 1, jpj 136 DO ji = 1, jpi 137 zcoefu = ( REAL(jpiglo+1)*0.5 - REAL(ji+nimpp-1) ) / ( REAL(jpiglo+1)*0.5 - 1. ) 138 zcoefv = ( REAL(jpjglo+1)*0.5 - REAL(jj+njmpp-1) ) / ( REAL(jpjglo+1)*0.5 - 1. ) 139 u_ice(ji,jj) = rn_uice * 1.5 * SIGN( 1., zcoefu ) * ABS( zcoefu ) * umask(ji,jj,1) 140 v_ice(ji,jj) = rn_vice * 1.5 * SIGN( 1., zcoefv ) * ABS( zcoefv ) * vmask(ji,jj,1) 141 END DO 142 END DO 137 143 ! --- 138 CALL ice_dyn_adv ( kt ) ! -- advection of ice144 CALL ice_dyn_adv ( kt ) ! -- advection of ice + correction on ice thickness 139 145 140 146 ! diagnostics: divergence at T points … … 150 156 DEALLOCATE( zdivu_i ) 151 157 158 CASE ( np_dynADV2D ) !== pure advection ==! (2D w prescribed velocities) 159 ALLOCATE( zdivu_i(jpi,jpj) ) 160 u_ice(:,:) = rn_uice * umask(:,:,1) 161 v_ice(:,:) = rn_vice * vmask(:,:,1) 162 !CALL RANDOM_NUMBER(u_ice(:,:)) ; u_ice(:,:) = u_ice(:,:) * 0.1 + rn_uice * 0.9 * umask(:,:,1) 163 !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1) 164 ! --- 165 CALL ice_dyn_adv ( kt ) ; CALL Hbig( zhi_max, zhs_max ) ! -- advection of ice + correction on ice thickness 166 167 ! diagnostics: divergence at T points 168 DO jj = 2, jpjm1 169 DO ji = 2, jpim1 170 zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 171 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj) 172 END DO 173 END DO 174 CALL lbc_lnk( zdivu_i, 'T', 1. ) 175 IF( iom_use('icediv') ) CALL iom_put( "icediv" , zdivu_i(:,:) ) 176 177 DEALLOCATE( zdivu_i ) 178 152 179 END SELECT 153 ! 154 IF( ln_timing ) CALL timing_stop('icedyn') 180 ! 181 ! controls 182 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icedyn', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 183 IF( ln_timing ) CALL timing_stop ('icedyn') ! timing 155 184 ! 156 185 END SUBROUTINE ice_dyn 157 186 158 187 159 SUBROUTINE Hbig( ph max )188 SUBROUTINE Hbig( phi_max, phs_max ) 160 189 !!------------------------------------------------------------------- 161 190 !! *** ROUTINE Hbig *** 162 191 !! 163 192 !! ** Purpose : Thickness correction in case advection scheme creates 164 !! abnormally tick ice 165 !! 166 !! ** Method : 1- check whether ice thickness resulting from advection is 167 !! larger than the surrounding 9-points before advection 168 !! and reduce it if a) divergence or b) convergence & at_i>0.8 169 !! 2- bound ice thickness with hi_max (99m) 193 !! abnormally tick ice or snow 194 !! 195 !! ** Method : 1- check whether ice thickness is larger than the surrounding 9-points 196 !! (before advection) and reduce it by adapting ice concentration 197 !! 2- check whether snow thickness is larger than the surrounding 9-points 198 !! (before advection) and reduce it by sending the excess in the ocean 199 !! 3- check whether snow load deplets the snow-ice interface below sea level$ 200 !! and reduce it by sending the excess in the ocean 201 !! 4- correct pond fraction to avoid a_ip > a_i 170 202 !! 171 203 !! ** input : Max thickness of the surrounding 9-points 172 204 !!------------------------------------------------------------------- 173 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: ph max ! max ice thick from surrounding 9-pts205 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi_max, phs_max ! max ice thick from surrounding 9-pts 174 206 ! 175 207 INTEGER :: ji, jj, jl ! dummy loop indices 176 REAL(wp) :: zh , zdv208 REAL(wp) :: zhi, zhs, zvs_excess, zfra 177 209 !!------------------------------------------------------------------- 178 210 ! … … 182 214 DO jj = 1, jpj 183 215 DO ji = 1, jpi 184 IF ( v_i(ji,jj,jl) > 0._wp ) THEN !-- bound to hmax216 IF ( v_i(ji,jj,jl) > 0._wp ) THEN 185 217 ! 186 zh = v_i (ji,jj,jl) / a_i(ji,jj,jl) 187 zdv = v_i(ji,jj,jl) - v_i_b(ji,jj,jl) 188 ! 189 IF ( ( zdv > 0.0 .AND. zh > phmax(ji,jj,jl) .AND. at_i_b(ji,jj) < 0.80 ) .OR. & 190 & ( zdv <= 0.0 .AND. zh > phmax(ji,jj,jl) ) ) THEN 191 a_i (ji,jj,jl) = v_i(ji,jj,jl) / MIN( phmax(ji,jj,jl), hi_max(jpl) ) !-- bound h_i to hi_max (99 m) 218 ! ! -- check h_i -- ! 219 ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i 220 zhi = v_i (ji,jj,jl) / a_i(ji,jj,jl) 221 !!clem zdv = v_i(ji,jj,jl) - v_i_b(ji,jj,jl) 222 !!clem IF ( ( zdv > 0.0 .AND. zh > phmax(ji,jj,jl) .AND. at_i_b(ji,jj) < 0.80 ) .OR. & 223 !!clem & ( zdv <= 0.0 .AND. zh > phmax(ji,jj,jl) ) ) THEN 224 IF( zhi > phi_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN 225 a_i(ji,jj,jl) = v_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) ) !-- bound h_i to hi_max (99 m) 192 226 ENDIF 193 227 ! 228 ! ! -- check h_s -- ! 229 ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean 230 zhs = v_s (ji,jj,jl) / a_i(ji,jj,jl) 231 IF( v_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN 232 zfra = a_i(ji,jj,jl) * phs_max(ji,jj,jl) / MAX( v_s(ji,jj,jl), epsi20 ) 233 ! 234 wfx_res(ji,jj) = wfx_res(ji,jj) + ( v_s(ji,jj,jl) - a_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * r1_rdtice 235 hfx_res(ji,jj) = hfx_res(ji,jj) - e_s(ji,jj,1,jl) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0 236 ! 237 e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * zfra 238 v_s(ji,jj,jl) = a_i(ji,jj,jl) * phs_max(ji,jj,jl) 239 ENDIF 240 ! 241 ! ! -- check snow load -- ! 242 ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean 243 ! this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin) 244 ! this imposed mini can artificially make the snow thickness very high (if concentration decreases drastically) 245 zvs_excess = MAX( 0._wp, v_s(ji,jj,jl) - v_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 246 IF( zvs_excess > 0._wp ) THEN 247 zfra = zvs_excess / MAX( v_s(ji,jj,jl), epsi20 ) 248 wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * r1_rdtice 249 hfx_res(ji,jj) = hfx_res(ji,jj) - e_s(ji,jj,1,jl) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0 250 ! 251 e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * zfra 252 v_s(ji,jj,jl) = v_s(ji,jj,jl) - zvs_excess 253 ENDIF 254 194 255 ENDIF 195 256 END DO … … 241 302 INTEGER :: ios, ioptio ! Local integer output status for namelist read 242 303 !! 243 NAMELIST/namdyn/ ln_dyn FULL, ln_dynRHGADV, ln_dynADV, rn_uice, rn_vice, &244 & rn_ishlat , &304 NAMELIST/namdyn/ ln_dynALL, ln_dynRHGADV, ln_dynADV1D, ln_dynADV2D, rn_uice, rn_vice, & 305 & rn_ishlat , & 245 306 & ln_landfast_L16, ln_landfast_home, rn_depfra, rn_icebfr, rn_lfrelax, rn_tensile 246 307 !!------------------------------------------------------------------- … … 259 320 WRITE(numout,*) '~~~~~~~~~~~~' 260 321 WRITE(numout,*) ' Namelist namdyn:' 261 WRITE(numout,*) ' Full ice dynamics (rhg + adv + ridge/raft + corr) ln_dyn FULL = ', ln_dynFULL322 WRITE(numout,*) ' Full ice dynamics (rhg + adv + ridge/raft + corr) ln_dynALL = ', ln_dynALL 262 323 WRITE(numout,*) ' No ridge/raft & No cor (rhg + adv) ln_dynRHGADV = ', ln_dynRHGADV 263 WRITE(numout,*) ' Advection only (rn_uvice + adv) ln_dynADV = ', ln_dynADV 324 WRITE(numout,*) ' Advection 1D only (Schar & Smolarkiewicz 1996) ln_dynADV1D = ', ln_dynADV1D 325 WRITE(numout,*) ' Advection 2D only (rn_uvice + adv) ln_dynADV2D = ', ln_dynADV2D 264 326 WRITE(numout,*) ' with prescribed velocity given by (u,v)_ice = (rn_uice,rn_vice) = (', rn_uice,',', rn_vice,')' 265 327 WRITE(numout,*) ' lateral boundary condition for sea ice dynamics rn_ishlat = ', rn_ishlat … … 275 337 ioptio = 0 276 338 ! !--- full dynamics (rheology + advection + ridging/rafting + correction) 277 IF( ln_dyn FULL ) THEN ; ioptio = ioptio + 1 ; nice_dyn = np_dynFULL; ENDIF339 IF( ln_dynALL ) THEN ; ioptio = ioptio + 1 ; nice_dyn = np_dynALL ; ENDIF 278 340 ! !--- dynamics without ridging/rafting and corr (rheology + advection) 279 341 IF( ln_dynRHGADV ) THEN ; ioptio = ioptio + 1 ; nice_dyn = np_dynRHGADV ; ENDIF 280 ! !--- advection only with prescribed ice velocities (from namelist) 281 IF( ln_dynADV ) THEN ; ioptio = ioptio + 1 ; nice_dyn = np_dynADV ; ENDIF 342 ! !--- advection 1D only - test case from Schar & Smolarkiewicz 1996 343 IF( ln_dynADV1D ) THEN ; ioptio = ioptio + 1 ; nice_dyn = np_dynADV1D ; ENDIF 344 ! !--- advection 2D only with prescribed ice velocities (from namelist) 345 IF( ln_dynADV2D ) THEN ; ioptio = ioptio + 1 ; nice_dyn = np_dynADV2D ; ENDIF 282 346 ! 283 347 IF( ioptio /= 1 ) CALL ctl_stop( 'ice_dyn_init: one and only one ice dynamics option has to be defined ' )
Note: See TracChangeset
for help on using the changeset viewer.
