Changeset 11081 for NEMO/branches/UKMO/NEMO_4.0_mirror/src/ICE
- Timestamp:
- 2019-06-06T16:11:54+02:00 (5 years ago)
- Location:
- NEMO/branches/UKMO/NEMO_4.0_mirror/src/ICE
- Files:
-
- 15 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0_mirror/src/ICE/icecor.F90
r10888 r11081 66 66 WRITE(numout,*) '~~~~~~~' 67 67 ENDIF 68 !69 68 ! !----------------------------------------------------- 70 ! ! ice thickness must exceed himin (for ice diff)!69 ! ! ice thickness must exceed himin (for temp. diff.) ! 71 70 ! !----------------------------------------------------- 72 71 WHERE( a_i(:,:,:) >= epsi20 ) ; h_i(:,:,:) = v_i(:,:,:) / a_i(:,:,:) … … 97 96 END DO 98 97 ENDIF 99 100 98 ! !----------------------------------------------------- 101 99 ! ! Rebin categories with thickness out of bounds ! -
NEMO/branches/UKMO/NEMO_4.0_mirror/src/ICE/icectl.F90
r10888 r11081 69 69 !! 70 70 REAL(wp) :: zv, zs, zt, zfs, zfv, zft 71 REAL(wp) :: zvmin, zamin, zamax 71 REAL(wp) :: zvmin, zamin, zamax, zeimin, zesmin, zsmin 72 72 REAL(wp) :: zvtrp, zetrp 73 73 REAL(wp) :: zarea, zv_sill, zs_sill, zt_sill … … 141 141 zetrp = glob_sum( 'icectl', ( diag_trp_ei + diag_trp_es ) * e1e2t ) * zconv 142 142 143 zvmin = glob_min( 'icectl', v_i ) 144 zamax = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 145 zamin = glob_min( 'icectl', a_i ) 143 zamax = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 144 zvmin = glob_min( 'icectl', v_i ) 145 zamin = glob_min( 'icectl', a_i ) 146 zsmin = glob_min( 'icectl', sv_i ) 147 zeimin = glob_min( 'icectl', SUM( e_i, dim=3 ) ) 148 zesmin = glob_min( 'icectl', SUM( e_s, dim=3 ) ) 146 149 147 150 ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice) … … 152 155 153 156 IF(lwp) THEN 154 IF ( ABS( zv ) > zv_sill ) WRITE(numout,*) 'violation volume [Mt/day] (',cd_routine,') = ',zv 155 IF ( ABS( zs ) > zs_sill ) WRITE(numout,*) 'violation saline [psu*Mt/day] (',cd_routine,') = ',zs 156 IF ( ABS( zt ) > zt_sill ) WRITE(numout,*) 'violation enthalpy [GW] (',cd_routine,') = ',zt 157 IF ( zvmin < -epsi10 ) WRITE(numout,*) 'violation v_i<0 [m] (',cd_routine,') = ',zvmin 158 IF ( zamax > MAX( rn_amax_n, rn_amax_s ) + epsi10 & 159 & .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' .AND. cd_routine /= 'Hbig' ) & 160 & WRITE(numout,*) 'violation a_i>amax (',cd_routine,') = ',zamax 161 IF ( zamin < -epsi10 ) WRITE(numout,*) 'violation a_i<0 (',cd_routine,') = ',zamin 162 !clem: the following check fails when using UMx advection scheme (see comments in icedyn_adv.F90) 157 ! check conservation issues 158 IF ( ABS( zv ) > zv_sill ) WRITE(numout,*) 'violation volume [Mt/day] (',cd_routine,') = ',zv 159 IF ( ABS( zs ) > zs_sill ) WRITE(numout,*) 'violation saline [psu*Mt/day] (',cd_routine,') = ',zs 160 IF ( ABS( zt ) > zt_sill ) WRITE(numout,*) 'violation enthalpy [GW] (',cd_routine,') = ',zt 161 ! check maximum ice concentration 162 IF ( zamax > MAX( rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' ) & 163 & WRITE(numout,*) 'violation a_i>amax (',cd_routine,') = ',zamax 164 ! check negative values 165 IF ( zvmin < 0. ) WRITE(numout,*) 'violation v_i<0 [m] (',cd_routine,') = ',zvmin 166 IF ( zamin < 0. ) WRITE(numout,*) 'violation a_i<0 (',cd_routine,') = ',zamin 167 IF ( zsmin < 0. ) WRITE(numout,*) 'violation s_i<0 (',cd_routine,') = ',zsmin 168 IF ( zeimin < 0. ) WRITE(numout,*) 'violation e_i<0 (',cd_routine,') = ',zeimin 169 IF ( zesmin < 0. ) WRITE(numout,*) 'violation e_s<0 (',cd_routine,') = ',zesmin 170 !clem: the following check fails (I think...) 163 171 ! IF ( ABS(zvtrp ) > zv_sill .AND. cd_routine == 'icedyn_adv' ) THEN 164 172 ! WRITE(numout,*) 'violation vtrp [Mt/day] (',cd_routine,') = ',zvtrp -
NEMO/branches/UKMO/NEMO_4.0_mirror/src/ICE/icedyn.F90
r10888 r11081 74 74 INTEGER, INTENT(in) :: kt ! ice time step 75 75 !! 76 INTEGER :: ji, jj , jl! dummy loop indices76 INTEGER :: ji, jj ! dummy loop indices 77 77 REAL(wp) :: zcoefu, zcoefv 78 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhi_max, zhs_max, zhip_max 79 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdivu_i 78 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdivu_i 80 79 !!-------------------------------------------------------------------- 81 80 ! … … 89 88 ENDIF 90 89 ! 91 IF( ln_landfast_home ) THEN !-- Landfast ice parameterization 92 tau_icebfr(:,:) = 0._wp 93 DO jl = 1, jpl 94 WHERE( h_i_b(:,:,jl) > ht_n(:,:) * rn_depfra ) tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 95 END DO 96 ENDIF 97 ! 98 ! !-- Record max of the surrounding 9-pts ice thick. (for CALL Hbig) 99 DO jl = 1, jpl 100 DO jj = 2, jpjm1 101 DO ji = fs_2, fs_jpim1 102 zhip_max(ji,jj,jl) = MAX( epsi20, h_ip_b(ji,jj,jl), h_ip_b(ji+1,jj ,jl), h_ip_b(ji ,jj+1,jl), & 103 & h_ip_b(ji-1,jj ,jl), h_ip_b(ji ,jj-1,jl), & 104 & h_ip_b(ji+1,jj+1,jl), h_ip_b(ji-1,jj-1,jl), & 105 & h_ip_b(ji+1,jj-1,jl), h_ip_b(ji-1,jj+1,jl) ) 106 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), & 107 & h_i_b (ji-1,jj ,jl), h_i_b (ji ,jj-1,jl), & 108 & h_i_b (ji+1,jj+1,jl), h_i_b (ji-1,jj-1,jl), & 109 & h_i_b (ji+1,jj-1,jl), h_i_b (ji-1,jj+1,jl) ) 110 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), & 111 & h_s_b (ji-1,jj ,jl), h_s_b (ji ,jj-1,jl), & 112 & h_s_b (ji+1,jj+1,jl), h_s_b (ji-1,jj-1,jl), & 113 & h_s_b (ji+1,jj-1,jl), h_s_b (ji-1,jj+1,jl) ) 114 END DO 115 END DO 116 END DO 117 CALL lbc_lnk_multi( 'icedyn', zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1., zhip_max(:,:,:), 'T', 1. ) 118 ! 119 ! 120 SELECT CASE( nice_dyn ) !-- Set which dynamics is running 90 ! retrieve thickness from volume for landfast param. and UMx advection scheme 91 WHERE( a_i(:,:,:) >= epsi20 ) 92 h_i(:,:,:) = v_i(:,:,:) / a_i_b(:,:,:) 93 h_s(:,:,:) = v_s(:,:,:) / a_i_b(:,:,:) 94 ELSEWHERE 95 h_i(:,:,:) = 0._wp 96 h_s(:,:,:) = 0._wp 97 END WHERE 98 ! 99 WHERE( a_ip(:,:,:) >= epsi20 ) 100 h_ip(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:) 101 ELSEWHERE 102 h_ip(:,:,:) = 0._wp 103 END WHERE 104 ! 105 ! 106 SELECT CASE( nice_dyn ) !-- Set which dynamics is running 121 107 122 108 CASE ( np_dynALL ) !== all dynamical processes ==! 123 CALL ice_dyn_rhg ( kt ) ! -- rheology 124 CALL ice_dyn_adv ( kt ) ; CALL Hbig( zhi_max, zhs_max, zhip_max ) ! -- advection of ice + correction on ice thickness 125 CALL ice_dyn_rdgrft( kt ) ! -- ridging/rafting 126 CALL ice_cor ( kt , 1 ) ! -- Corrections 127 109 ! 110 CALL ice_dyn_rhg ( kt ) ! -- rheology 111 CALL ice_dyn_adv ( kt ) ! -- advection of ice 112 CALL ice_dyn_rdgrft( kt ) ! -- ridging/rafting 113 CALL ice_cor ( kt , 1 ) ! -- Corrections 114 ! 128 115 CASE ( np_dynRHGADV ) !== no ridge/raft & no corrections ==! 129 CALL ice_dyn_rhg ( kt ) ! -- rheology 130 CALL ice_dyn_adv ( kt ) ; CALL Hbig( zhi_max, zhs_max, zhip_max ) ! -- advection of ice + correction on ice thickness 131 CALL Hpiling ! -- simple pile-up (replaces ridging/rafting) 132 116 ! 117 CALL ice_dyn_rhg ( kt ) ! -- rheology 118 CALL ice_dyn_adv ( kt ) ! -- advection of ice 119 CALL Hpiling ! -- simple pile-up (replaces ridging/rafting) 120 CALL ice_var_zapsmall ! -- zap small areas 121 ! 133 122 CASE ( np_dynADV1D ) !== pure advection ==! (1D) 134 ALLOCATE( zdivu_i(jpi,jpj) )123 ! 135 124 ! --- monotonicity test from Schar & Smolarkiewicz 1996 --- ! 136 125 ! CFL = 0.5 at a distance from the bound of 1/6 of the basin length … … 145 134 END DO 146 135 ! --- 147 CALL ice_dyn_adv ( kt ) ! -- advection of ice 148 149 ! diagnostics: divergence at T points 150 DO jj = 2, jpjm1 151 DO ji = 2, jpim1 152 zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 153 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj) 154 END DO 155 END DO 156 CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. ) 157 IF( iom_use('icediv') ) CALL iom_put( "icediv" , zdivu_i(:,:) ) 158 159 DEALLOCATE( zdivu_i ) 160 136 CALL ice_dyn_adv ( kt ) ! -- advection of ice 137 ! 161 138 CASE ( np_dynADV2D ) !== pure advection ==! (2D w prescribed velocities) 162 ALLOCATE( zdivu_i(jpi,jpj) )139 ! 163 140 u_ice(:,:) = rn_uice * umask(:,:,1) 164 141 v_ice(:,:) = rn_vice * vmask(:,:,1) … … 166 143 !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1) 167 144 ! --- 168 CALL ice_dyn_adv ( kt ) ! -- advection of ice 169 170 ! diagnostics: divergence at T points 171 DO jj = 2, jpjm1 172 DO ji = 2, jpim1 173 zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 174 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj) 145 CALL ice_dyn_adv ( kt ) ! -- advection of ice 146 147 END SELECT 148 ! 149 ! 150 ! diagnostics: divergence at T points 151 IF( iom_use('icediv') ) THEN 152 ! 153 SELECT CASE( nice_dyn ) 154 155 CASE ( np_dynADV1D , np_dynADV2D ) 156 157 ALLOCATE( zdivu_i(jpi,jpj) ) 158 DO jj = 2, jpjm1 159 DO ji = 2, jpim1 160 zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 161 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj) 162 END DO 175 163 END DO 176 END DO177 CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1.)178 IF( iom_use('icediv') ) CALL iom_put( "icediv" , zdivu_i(:,:))179 180 DEALLOCATE( zdivu_i )181 182 END SELECT164 CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. ) 165 CALL iom_put( "icediv" , zdivu_i(:,:) ) 166 DEALLOCATE( zdivu_i ) 167 168 END SELECT 169 ! 170 ENDIF 183 171 ! 184 172 ! controls … … 188 176 189 177 190 SUBROUTINE Hbig( phi_max, phs_max, phip_max )191 !!-------------------------------------------------------------------192 !! *** ROUTINE Hbig ***193 !!194 !! ** Purpose : Thickness correction in case advection scheme creates195 !! abnormally tick ice or snow196 !!197 !! ** Method : 1- check whether ice thickness is larger than the surrounding 9-points198 !! (before advection) and reduce it by adapting ice concentration199 !! 2- check whether snow thickness is larger than the surrounding 9-points200 !! (before advection) and reduce it by sending the excess in the ocean201 !! 3- check whether snow load deplets the snow-ice interface below sea level$202 !! and reduce it by sending the excess in the ocean203 !! 4- correct pond fraction to avoid a_ip > a_i204 !!205 !! ** input : Max thickness of the surrounding 9-points206 !!-------------------------------------------------------------------207 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi_max, phs_max, phip_max ! max ice thick from surrounding 9-pts208 !209 INTEGER :: ji, jj, jl ! dummy loop indices210 REAL(wp) :: zhip, zhi, zhs, zvs_excess, zfra211 !!-------------------------------------------------------------------212 ! controls213 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'Hbig', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation214 !215 CALL ice_var_zapsmall !-- zap small areas216 !217 DO jl = 1, jpl218 DO jj = 1, jpj219 DO ji = 1, jpi220 IF ( v_i(ji,jj,jl) > 0._wp ) THEN221 !222 ! ! -- check h_ip -- !223 ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip224 IF( ln_pnd_H12 .AND. v_ip(ji,jj,jl) > 0._wp ) THEN225 zhip = v_ip(ji,jj,jl) / MAX( epsi20, a_ip(ji,jj,jl) )226 IF( zhip > phip_max(ji,jj,jl) .AND. a_ip(ji,jj,jl) < 0.15 ) THEN227 a_ip(ji,jj,jl) = v_ip(ji,jj,jl) / phip_max(ji,jj,jl)228 ENDIF229 ENDIF230 !231 ! ! -- check h_i -- !232 ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i233 zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)234 IF( zhi > phi_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN235 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)236 ENDIF237 !238 ! ! -- check h_s -- !239 ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean240 zhs = v_s(ji,jj,jl) / a_i(ji,jj,jl)241 IF( v_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN242 zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 )243 !244 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_rdtice245 hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0246 !247 e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra248 v_s(ji,jj,jl) = a_i(ji,jj,jl) * phs_max(ji,jj,jl)249 ENDIF250 !251 ! ! -- check snow load -- !252 ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean253 ! this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin)254 ! this imposed mini can artificially make the snow very thick (if concentration decreases drastically)255 zvs_excess = MAX( 0._wp, v_s(ji,jj,jl) - v_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos )256 IF( zvs_excess > 0._wp ) THEN257 zfra = ( v_s(ji,jj,jl) - zvs_excess ) / MAX( v_s(ji,jj,jl), epsi20 )258 wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * r1_rdtice259 hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0260 !261 e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra262 v_s(ji,jj,jl) = v_s(ji,jj,jl) - zvs_excess263 ENDIF264 265 ENDIF266 END DO267 END DO268 END DO269 ! !-- correct pond fraction to avoid a_ip > a_i270 WHERE( a_ip(:,:,:) > a_i(:,:,:) ) a_ip(:,:,:) = a_i(:,:,:)271 !272 ! controls273 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'Hbig', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation274 !275 END SUBROUTINE Hbig276 277 278 178 SUBROUTINE Hpiling 279 179 !!------------------------------------------------------------------- … … 290 190 ! controls 291 191 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'Hpiling', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 292 !293 CALL ice_var_zapsmall !-- zap small areas294 192 ! 295 193 at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) … … 337 235 WRITE(numout,*) '~~~~~~~~~~~~' 338 236 WRITE(numout,*) ' Namelist namdyn:' 339 WRITE(numout,*) ' Full ice dynamics (rhg + adv + ridge/raft + corr) 340 WRITE(numout,*) ' No ridge/raft & No cor (rhg + adv) 341 WRITE(numout,*) ' Advection 1D only (Schar & Smolarkiewicz 1996) 342 WRITE(numout,*) ' Advection 2D only (rn_uvice + adv) 343 WRITE(numout,*) ' with prescribed velocity given by (u,v)_ice = (rn_uice,rn_vice) = (', rn_uice,',',rn_vice,')'344 WRITE(numout,*) ' lateral boundary condition for sea ice dynamics 345 WRITE(numout,*) ' Landfast: param from Lemieux 2016 346 WRITE(numout,*) ' Landfast: param from home made 347 WRITE(numout,*) ' fraction of ocean depth that ice must reach 348 WRITE(numout,*) ' maximum bottom stress per unit area of contact 349 WRITE(numout,*) ' relax time scale (s-1) to reach static friction 350 WRITE(numout,*) ' isotropic tensile strength 237 WRITE(numout,*) ' Full ice dynamics (rhg + adv + ridge/raft + corr) ln_dynALL = ', ln_dynALL 238 WRITE(numout,*) ' No ridge/raft & No cor (rhg + adv) ln_dynRHGADV = ', ln_dynRHGADV 239 WRITE(numout,*) ' Advection 1D only (Schar & Smolarkiewicz 1996) ln_dynADV1D = ', ln_dynADV1D 240 WRITE(numout,*) ' Advection 2D only (rn_uvice + adv) ln_dynADV2D = ', ln_dynADV2D 241 WRITE(numout,*) ' with prescribed velocity given by (u,v)_ice = (rn_uice,rn_vice) = (', rn_uice,',',rn_vice,')' 242 WRITE(numout,*) ' lateral boundary condition for sea ice dynamics rn_ishlat = ', rn_ishlat 243 WRITE(numout,*) ' Landfast: param from Lemieux 2016 ln_landfast_L16 = ', ln_landfast_L16 244 WRITE(numout,*) ' Landfast: param from home made ln_landfast_home= ', ln_landfast_home 245 WRITE(numout,*) ' fraction of ocean depth that ice must reach rn_depfra = ', rn_depfra 246 WRITE(numout,*) ' maximum bottom stress per unit area of contact rn_icebfr = ', rn_icebfr 247 WRITE(numout,*) ' relax time scale (s-1) to reach static friction rn_lfrelax = ', rn_lfrelax 248 WRITE(numout,*) ' isotropic tensile strength rn_tensile = ', rn_tensile 351 249 WRITE(numout,*) 352 250 ENDIF -
NEMO/branches/UKMO/NEMO_4.0_mirror/src/ICE/icedyn_adv.F90
r10888 r11081 64 64 !!---------------------------------------------------------------------- 65 65 INTEGER, INTENT(in) :: kt ! number of iteration 66 !67 INTEGER :: jl ! dummy loop indice68 REAL(wp), DIMENSION(jpi,jpj) :: zmask ! fraction of time step with v_i < 069 66 !!--------------------------------------------------------------------- 70 67 ! 71 IF( ln_timing ) CALL timing_start('icedyn_adv') 68 ! controls 69 IF( ln_timing ) CALL timing_start('icedyn_adv') ! timing 70 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icedyn_adv', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 72 71 ! 73 72 IF( kt == nit000 .AND. lwp ) THEN … … 76 75 WRITE(numout,*) '~~~~~~~~~~~' 77 76 ENDIF 78 79 ! conservation test 80 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icedyn_adv', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 81 82 !---------- 83 ! Advection 84 !---------- 77 ! 78 !---------------! 79 !== Advection ==! 80 !---------------! 85 81 SELECT CASE( nice_adv ) 86 82 ! !-----------------------! 87 83 CASE( np_advUMx ) ! ULTIMATE-MACHO scheme ! 88 84 ! !-----------------------! 89 CALL ice_dyn_adv_umx( nn_UMx, kt, u_ice, v_ice, ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 90 ! !-----------------------! 85 CALL ice_dyn_adv_umx( nn_UMx, kt, u_ice, v_ice, h_i, h_s, h_ip, & 86 & ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 87 ! !-----------------------! 91 88 CASE( np_advPRA ) ! PRATHER scheme ! 92 89 ! !-----------------------! 93 CALL ice_dyn_adv_pra( kt, u_ice, v_ice, ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 90 CALL ice_dyn_adv_pra( kt, u_ice, v_ice, & 91 & ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 94 92 END SELECT 95 96 !----------------------------97 ! Debug the advection schemes98 !----------------------------99 ! clem: At least one advection scheme above is not strictly positive => UMx100 ! In Prather, I am not sure if the fields are bounded by 0 or not (it seems yes)101 ! In UMx , advected fields are not perfectly bounded and negative values can appear.102 ! These values are usually very small but in some occasions they can also be non-negligible103 ! Therefore one needs to bound the advected fields by 0 (though this is not a clean fix)104 !105 ! record the negative values resulting from UMx106 zmask(:,:) = 0._wp ! keep the init to 0 here107 DO jl = 1, jpl108 WHERE( v_i(:,:,jl) < 0._wp ) zmask(:,:) = 1._wp109 END DO110 IF( iom_use('iceneg_pres') ) CALL iom_put("iceneg_pres", zmask ) ! fraction of time step with v_i < 0111 IF( iom_use('iceneg_volu') ) CALL iom_put("iceneg_volu", SUM(MIN( v_i, 0. ), dim=3 ) ) ! negative ice volume (only)112 IF( iom_use('iceneg_hfx' ) ) CALL iom_put("iceneg_hfx" , ( SUM(SUM( MIN( e_i(:,:,1:nlay_i,:), 0. ) & ! negative ice heat content (only)113 & , dim=4 ), dim=3 ) )* r1_rdtice ) ! -- eq. heat flux [W/m2]114 !115 ! ==> conservation is ensured by calling this routine below,116 ! however the global ice volume is then changed by advection (but errors are small)117 CALL ice_var_zapneg( ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )118 93 119 94 !------------ -
NEMO/branches/UKMO/NEMO_4.0_mirror/src/ICE/icedyn_adv_umx.F90
r10888 r11081 11 11 !! 'key_si3' SI3 sea-ice model 12 12 !!---------------------------------------------------------------------- 13 !! ice_dyn_adv_umx : update the tracer trend with the 3D advection trends using a TVD scheme13 !! ice_dyn_adv_umx : update the tracer fields 14 14 !! ultimate_x(_y) : compute a tracer value at velocity points using ULTIMATE scheme at various orders 15 !! macho : ???16 !! nonosc_ice : compute monotonic tracer fluxes bya non-oscillatory algorithm15 !! macho : compute the fluxes 16 !! nonosc_ice : limit the fluxes using a non-oscillatory algorithm 17 17 !!---------------------------------------------------------------------- 18 18 USE phycst ! physical constant … … 32 32 33 33 PUBLIC ice_dyn_adv_umx ! called by icedyn_adv.F90 34 35 REAL(wp) :: z1_6 = 1._wp / 6._wp ! =1/6 36 REAL(wp) :: z1_120 = 1._wp / 120._wp ! =1/120 37 38 ! limiter: 1=nonosc_ice, 2=superbee, 3=h3(rachid) 39 INTEGER :: kn_limiter = 1 40 41 ! if T interpolated at u/v points is negative, then interpolate T at u/v points using the upstream scheme 42 ! clem: if set to true, the 2D test case "diagonal advection" does not work (I do not understand why) 43 ! but in realistic cases, it avoids having very negative ice temperature (-50) at low ice concentration 44 LOGICAL :: ll_neg = .TRUE. 45 46 ! alternate directions for upstream 47 LOGICAL :: ll_upsxy = .TRUE. 48 49 ! alternate directions for high order 50 LOGICAL :: ll_hoxy = .TRUE. 51 52 ! prelimiter: use it to avoid overshoot in H 53 ! clem: if set to true, the 2D test case "diagnoal advection" does not work (I do not understand why) 54 LOGICAL :: ll_prelimiter_zalesak = .FALSE. ! from: Zalesak(1979) eq. 14 => better for 1D. Not well defined in 2D 55 56 34 ! 35 INTEGER, PARAMETER :: np_advS = 1 ! advection for S and T: dVS/dt = -div( uVS ) => np_advS = 1 36 ! or dVS/dt = -div( uA * uHS / u ) => np_advS = 2 37 ! or dVS/dt = -div( uV * uS / u ) => np_advS = 3 38 INTEGER, PARAMETER :: np_limiter = 1 ! limiter: 1 = nonosc 39 ! 2 = superbee 40 ! 3 = h3 41 LOGICAL :: ll_upsxy = .TRUE. ! alternate directions for upstream 42 LOGICAL :: ll_hoxy = .TRUE. ! alternate directions for high order 43 LOGICAL :: ll_neg = .TRUE. ! if T interpolated at u/v points is negative or v_i < 1.e-6 44 ! then interpolate T at u/v points using the upstream scheme 45 LOGICAL :: ll_prelim = .FALSE. ! prelimiter from: Zalesak(1979) eq. 14 => not well defined in 2D 46 ! 47 REAL(wp) :: z1_6 = 1._wp / 6._wp ! =1/6 48 REAL(wp) :: z1_120 = 1._wp / 120._wp ! =1/120 49 ! 50 INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: imsk_small, jmsk_small 51 ! 57 52 !! * Substitutions 58 53 # include "vectopt_loop_substitute.h90" … … 64 59 CONTAINS 65 60 66 SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, &61 SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip, & 67 62 & pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 68 63 !!---------------------------------------------------------------------- … … 79 74 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pu_ice ! ice i-velocity 80 75 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pv_ice ! ice j-velocity 76 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_i ! ice thickness 77 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_s ! snw thickness 78 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_ip ! ice pond thickness 81 79 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pato_i ! open water area 82 80 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i ! ice volume … … 93 91 INTEGER :: icycle ! number of sub-timestep for the advection 94 92 REAL(wp) :: zamsk ! 1 if advection of concentration, 0 if advection of other tracers 95 REAL(wp) :: zdt 96 REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! send zcflnow and receive zcflprv97 REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx, zcu_box, zcv_box 93 REAL(wp) :: zdt, zvi_cen 94 REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! for global communication 95 REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx, zcu_box, zcv_box 98 96 REAL(wp), DIMENSION(jpi,jpj) :: zati1, zati2 99 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zua_ho, zva_ho 100 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_ai, z1_aip 101 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhvar 97 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zu_cat, zv_cat 98 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zua_ho, zva_ho, zua_ups, zva_ups 99 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_ai , z1_aip, zhvar 100 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhi_max, zhs_max, zhip_max 101 ! 102 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs 102 103 !!---------------------------------------------------------------------- 103 104 ! 104 105 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme' 105 106 ! 106 ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- ! 107 ! When needed, the advection split is applied at the next time-step in order to avoid blocking global comm. 108 ! ...this should not affect too much the stability... Was ok on the tests we did... 107 ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- ! 108 DO jl = 1, jpl 109 DO jj = 2, jpjm1 110 DO ji = fs_2, fs_jpim1 111 zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj ,jl), ph_ip(ji ,jj+1,jl), & 112 & ph_ip(ji-1,jj ,jl), ph_ip(ji ,jj-1,jl), & 113 & ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 114 & ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 115 zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj ,jl), ph_i (ji ,jj+1,jl), & 116 & ph_i (ji-1,jj ,jl), ph_i (ji ,jj-1,jl), & 117 & ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 118 & ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 119 zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj ,jl), ph_s (ji ,jj+1,jl), & 120 & ph_s (ji-1,jj ,jl), ph_s (ji ,jj-1,jl), & 121 & ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 122 & ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 123 END DO 124 END DO 125 END DO 126 CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1. ) 127 ! 128 ! 129 ! --- If ice drift is too fast, use subtime steps for advection (CFL test for stability) --- ! 130 ! Note: the advection split is applied at the next time-step in order to avoid blocking global comm. 131 ! this should not affect too much the stability 109 132 zcflnow(1) = MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) 110 133 zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) … … 116 139 ELSE ; icycle = 1 117 140 ENDIF 118 119 141 zdt = rdt_ice / REAL(icycle) 120 142 … … 122 144 zudy(:,:) = pu_ice(:,:) * e2u(:,:) 123 145 zvdx(:,:) = pv_ice(:,:) * e1v(:,:) 124 146 ! 147 ! setup transport for each ice cat 148 DO jl = 1, jpl 149 zu_cat(:,:,jl) = zudy(:,:) 150 zv_cat(:,:,jl) = zvdx(:,:) 151 END DO 152 ! 125 153 ! --- define velocity for advection: u*grad(H) --- ! 126 154 DO jj = 2, jpjm1 … … 154 182 END WHERE 155 183 ! 156 ! set u*a=u for advection of A only 157 DO jl = 1, jpl 158 zua_ho(:,:,jl) = zudy(:,:) 159 zva_ho(:,:,jl) = zvdx(:,:) 160 END DO 161 184 ! setup a mask where advection will be upstream 185 IF( ll_neg ) THEN 186 IF( .NOT. ALLOCATED(imsk_small) ) ALLOCATE( imsk_small(jpi,jpj,jpl) ) 187 IF( .NOT. ALLOCATED(jmsk_small) ) ALLOCATE( jmsk_small(jpi,jpj,jpl) ) 188 DO jl = 1, jpl 189 DO jj = 1, jpjm1 190 DO ji = 1, jpim1 191 zvi_cen = 0.5_wp * ( pv_i(ji+1,jj,jl) + pv_i(ji,jj,jl) ) 192 IF( zvi_cen < epsi06) THEN ; imsk_small(ji,jj,jl) = 0 193 ELSE ; imsk_small(ji,jj,jl) = 1 ; ENDIF 194 zvi_cen = 0.5_wp * ( pv_i(ji,jj+1,jl) + pv_i(ji,jj,jl) ) 195 IF( zvi_cen < epsi06) THEN ; jmsk_small(ji,jj,jl) = 0 196 ELSE ; jmsk_small(ji,jj,jl) = 1 ; ENDIF 197 END DO 198 END DO 199 END DO 200 ENDIF 201 ! 202 ! ----------------------- ! 203 ! ==> start advection <== ! 204 ! ----------------------- ! 205 ! 206 !== Ice area ==! 162 207 zamsk = 1._wp 163 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_i, pa_i, zua_ho, zva_ho ) !-- Ice area 164 zamsk = 0._wp 165 ! 166 zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:) 167 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_i ) !-- Ice volume 168 ! 169 zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:) 170 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_s ) !-- Snw volume 171 ! 172 zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:) 173 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, psv_i ) !-- Salt content 174 ! 175 DO jk = 1, nlay_i 176 zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:) 177 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_i(:,:,jk,:) ) !-- Ice heat content 178 END DO 179 ! 180 DO jk = 1, nlay_s 181 zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:) 182 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_s(:,:,jk,:) ) !-- Snw heat content 183 END DO 184 ! 185 IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN !-- Ice Age 186 ! clem: in theory we should use the formulation below to advect the ice age, but the code is unable to deal with 187 ! fields that do not depend on volume (here oa_i depends on concentration). It creates abnormal ages that 188 ! spread into the domain. Therefore we cheat and consider that ice age should be advected as ice concentration 189 !!zhvar(:,:,:) = poa_i(:,:,:) * z1_ai(:,:,:) 190 !!CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, poa_i ) 191 ! set u*a=u for advection of ice age 192 DO jl = 1, jpl 193 zua_ho(:,:,jl) = zudy(:,:) 194 zva_ho(:,:,jl) = zvdx(:,:) 195 END DO 208 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zu_cat , zv_cat , zcu_box, zcv_box, & 209 & pa_i, pa_i, zua_ups, zva_ups, zua_ho , zva_ho ) 210 ! 211 ! ! --------------------------------- ! 212 IF( np_advS == 1 ) THEN ! -- advection form: -div( uVS ) -- ! 213 ! ! --------------------------------- ! 214 zamsk = 0._wp 215 !== Ice volume ==! 216 zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:) 217 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 218 & zhvar, pv_i, zua_ups, zva_ups ) 219 !== Snw volume ==! 220 zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:) 221 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 222 & zhvar, pv_s, zua_ups, zva_ups ) 223 ! 196 224 zamsk = 1._wp 197 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, poa_i, poa_i ) 225 !== Salt content ==! 226 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, & 227 & psv_i, psv_i ) 228 !== Ice heat content ==! 229 DO jk = 1, nlay_i 230 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, & 231 & pe_i(:,:,jk,:), pe_i(:,:,jk,:) ) 232 END DO 233 !== Snw heat content ==! 234 DO jk = 1, nlay_s 235 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, & 236 & pe_s(:,:,jk,:), pe_s(:,:,jk,:) ) 237 END DO 238 ! 239 ! ! ------------------------------------------ ! 240 ELSEIF( np_advS == 2 ) THEN ! -- advection form: -div( uA * uHS / u ) -- ! 241 ! ! ------------------------------------------ ! 198 242 zamsk = 0._wp 243 !== Ice volume ==! 244 zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:) 245 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 246 & zhvar, pv_i, zua_ups, zva_ups ) 247 !== Snw volume ==! 248 zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:) 249 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 250 & zhvar, pv_s, zua_ups, zva_ups ) 251 !== Salt content ==! 252 zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:) 253 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 254 & zhvar, psv_i, zua_ups, zva_ups ) 255 !== Ice heat content ==! 256 DO jk = 1, nlay_i 257 zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:) 258 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho, zva_ho, zcu_box, zcv_box, & 259 & zhvar, pe_i(:,:,jk,:), zua_ups, zva_ups ) 260 END DO 261 !== Snw heat content ==! 262 DO jk = 1, nlay_s 263 zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:) 264 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho, zva_ho, zcu_box, zcv_box, & 265 & zhvar, pe_s(:,:,jk,:), zua_ups, zva_ups ) 266 END DO 267 ! 268 ! ! ----------------------------------------- ! 269 ELSEIF( np_advS == 3 ) THEN ! -- advection form: -div( uV * uS / u ) -- ! 270 ! ! ----------------------------------------- ! 271 zamsk = 0._wp 272 ! 273 ALLOCATE( zuv_ho (jpi,jpj,jpl), zvv_ho (jpi,jpj,jpl), & 274 & zuv_ups(jpi,jpj,jpl), zvv_ups(jpi,jpj,jpl), z1_vi(jpi,jpj,jpl), z1_vs(jpi,jpj,jpl) ) 275 ! 276 ! inverse of Vi 277 WHERE( pv_i(:,:,:) >= epsi20 ) ; z1_vi(:,:,:) = 1._wp / pv_i(:,:,:) 278 ELSEWHERE ; z1_vi(:,:,:) = 0. 279 END WHERE 280 ! inverse of Vs 281 WHERE( pv_s(:,:,:) >= epsi20 ) ; z1_vs(:,:,:) = 1._wp / pv_s(:,:,:) 282 ELSEWHERE ; z1_vs(:,:,:) = 0. 283 END WHERE 284 ! 285 ! It is important to first calculate the ice fields and then the snow fields (because we use the same arrays) 286 ! 287 !== Ice volume ==! 288 zuv_ups = zua_ups 289 zvv_ups = zva_ups 290 zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:) 291 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 292 & zhvar, pv_i, zuv_ups, zvv_ups, zuv_ho , zvv_ho ) 293 !== Salt content ==! 294 zhvar(:,:,:) = psv_i(:,:,:) * z1_vi(:,:,:) 295 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zuv_ho , zvv_ho , zcu_box, zcv_box, & 296 & zhvar, psv_i, zuv_ups, zvv_ups ) 297 !== Ice heat content ==! 298 DO jk = 1, nlay_i 299 zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_vi(:,:,:) 300 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zuv_ho, zvv_ho, zcu_box, zcv_box, & 301 & zhvar, pe_i(:,:,jk,:), zuv_ups, zvv_ups ) 302 END DO 303 !== Snow volume ==! 304 zuv_ups = zua_ups 305 zvv_ups = zva_ups 306 zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:) 307 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 308 & zhvar, pv_s, zuv_ups, zvv_ups, zuv_ho , zvv_ho ) 309 !== Snw heat content ==! 310 DO jk = 1, nlay_s 311 zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_vs(:,:,:) 312 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zuv_ho, zvv_ho, zcu_box, zcv_box, & 313 & zhvar, pe_s(:,:,jk,:), zuv_ups, zvv_ups ) 314 END DO 315 ! 316 DEALLOCATE( zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs ) 317 ! 199 318 ENDIF 200 319 ! 201 IF ( ln_pnd_H12 ) THEN !-- melt ponds 202 ! set u*a=u for advection of Ap only 203 DO jl = 1, jpl 204 zua_ho(:,:,jl) = zudy(:,:) 205 zva_ho(:,:,jl) = zvdx(:,:) 206 END DO 207 ! 320 !== Ice age ==! 321 IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN 208 322 zamsk = 1._wp 209 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_ip, pa_ip, zua_ho, zva_ho ) ! fraction 323 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, & 324 & poa_i, poa_i ) 325 ENDIF 326 ! 327 !== melt ponds ==! 328 IF ( ln_pnd_H12 ) THEN 329 ! fraction 330 zamsk = 1._wp 331 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat , zv_cat , zcu_box, zcv_box, & 332 & pa_ip, pa_ip, zua_ups, zva_ups, zua_ho , zva_ho ) 333 ! volume 210 334 zamsk = 0._wp 211 !212 335 zhvar(:,:,:) = pv_ip(:,:,:) * z1_aip(:,:,:) 213 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_ip ) ! volume 336 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 337 & zhvar, pv_ip, zua_ups, zva_ups ) 214 338 ENDIF 215 339 ! 340 !== Open water area ==! 216 341 zati2(:,:) = SUM( pa_i(:,:,:), dim=3 ) 217 342 DO jj = 2, jpjm1 218 343 DO ji = fs_2, fs_jpim1 219 pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & !-- Open water area344 pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & 220 345 & - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 221 346 END DO 222 347 END DO 223 CALL lbc_lnk( 'icedyn_adv_umx', pato_i(:,:), 'T', 1. ) 224 ! 348 CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T', 1. ) 349 ! 350 ! 351 ! --- Ensure non-negative fields and in-bound thicknesses --- ! 352 ! Remove negative values (conservation is ensured) 353 ! (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 354 CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 355 ! 356 ! Make sure ice thickness is not too big 357 ! (because ice thickness can be too large where ice concentration is very small) 358 CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 359 225 360 END DO 226 361 ! … … 228 363 229 364 230 SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pua_ho, pva_ho ) 365 SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, & 366 & pt, ptc, pua_ups, pva_ups, pua_ho, pva_ho ) 231 367 !!---------------------------------------------------------------------- 232 368 !! *** ROUTINE adv_umx *** … … 235 371 !! tracers and add it to the general trend of tracer equations 236 372 !! 237 !! ** Method : - calculate upstream fluxes and upstream solution for tracer H373 !! ** Method : - calculate upstream fluxes and upstream solution for tracers V/A(=H) etc 238 374 !! - calculate tracer H at u and v points (Ultimate) 239 !! - calculate the high order fluxes using alterning directions (Macho ?)375 !! - calculate the high order fluxes using alterning directions (Macho) 240 376 !! - apply a limiter on the fluxes (nonosc_ice) 241 !! - convert this tracer flux to a tracer content flux (uH -> uV) 242 !! - calculate the high order solution for tracer content V 377 !! - convert this tracer flux to a "volume" flux (uH -> uV) 378 !! - apply a limiter a second time on the volumes fluxes (nonosc_ice) 379 !! - calculate the high order solution for V 243 380 !! 244 !! ** Action : solve 2 equations => a) da/dt = -div(ua) 245 !! b) dV/dt = -div(uV) using dH/dt = -u.grad(H) 246 !! in eq. b), - fluxes uH are evaluated (with UMx) and limited (with nonosc_ice). This step is necessary to get a good H. 247 !! - then we convert this flux to a "volume" flux this way => uH*ua/u 248 !! where ua is the flux from eq. a) 249 !! - at last we estimate dV/dt = -div(uH*ua/u) 381 !! ** Action : solve 3 equations => a) dA/dt = -div(uA) 382 !! b) dV/dt = -div(uV) using dH/dt = -u.grad(H) 383 !! c) dVS/dt = -div(uVS) using either dHS/dt = -u.grad(HS) or dS/dt = -u.grad(S) 250 384 !! 251 !! ** Note : - this method can lead to small negative V (since we only limit H) => corrected in icedyn_adv.F90 conserving mass etc. 252 !! - negative tracers at u-v points can also occur from the Ultimate scheme (usually at the ice edge) and the solution for now 253 !! is to apply an upstream scheme when it occurs. A better solution would be to degrade the order of 254 !! the scheme automatically by applying a mask of the ice cover inside Ultimate (not done). 385 !! in eq. b), - fluxes uH are evaluated (with UMx) and limited with nonosc_ice. This step is necessary to get a good H. 386 !! - then we convert this flux to a "volume" flux this way => uH * uA / u 387 !! where uA is the flux from eq. a) 388 !! this "volume" flux is also limited with nonosc_ice (otherwise overshoots can occur) 389 !! - at last we estimate dV/dt = -div(uH * uA / u) 390 !! 391 !! in eq. c), one can solve the equation for S (ln_advS=T), then dVS/dt = -div(uV * uS / u) 392 !! or for HS (ln_advS=F), then dVS/dt = -div(uA * uHS / u) 393 !! 394 !! ** Note : - this method can lead to tiny negative V (-1.e-20) => set it to 0 while conserving mass etc. 395 !! - At the ice edge, Ultimate scheme can lead to: 396 !! 1) negative interpolated tracers at u-v points 397 !! 2) non-zero interpolated tracers at u-v points eventhough there is no ice and velocity is outward 398 !! Solution for 1): apply an upstream scheme when it occurs. A better solution would be to degrade the order of 399 !! the scheme automatically by applying a mask of the ice cover inside Ultimate (not done). 400 !! Solution for 2): we set it to 0 in this case 255 401 !! - Eventhough 1D tests give very good results (typically the one from Schar & Smolarkiewiecz), the 2D is less good. 256 402 !! Large values of H can appear for very small ice concentration, and when it does it messes the things up since we 257 !! work on H (and not V). It probably comes from the prelimiter of zalesak which is coded for 1D and not 2D.403 !! work on H (and not V). It is partly related to the multi-category approach 258 404 !! Therefore, after advection we limit the thickness to the largest value of the 9-points around (only if ice 259 !! concentration is small). 260 !! To-do: expand the prelimiter from zalesak to make it work in 2D 261 !!---------------------------------------------------------------------- 262 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) 263 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 264 INTEGER , INTENT(in ) :: jt ! number of sub-iteration 265 INTEGER , INTENT(in ) :: kt ! number of iteration 266 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 267 REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu , pv ! 2 ice velocity components => u*e2 268 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: puc , pvc ! 2 ice velocity components => u*e2 or u*a*e2u 269 REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pubox, pvbox ! upstream velocity 270 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pt ! tracer field 271 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: ptc ! tracer content field 272 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out), OPTIONAL :: pua_ho, pva_ho ! high order u*a fluxes 405 !! concentration is small). Since we do not limit S and T, large values can occur at the edge but it does not really matter 406 !! since sv_i and e_i are still good. 407 !!---------------------------------------------------------------------- 408 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) 409 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 410 INTEGER , INTENT(in ) :: jt ! number of sub-iteration 411 INTEGER , INTENT(in ) :: kt ! number of iteration 412 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 413 REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu , pv ! 2 ice velocity components => u*e2 414 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: puc , pvc ! 2 ice velocity components => u*e2 or u*a*e2u 415 REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pubox, pvbox ! upstream velocity 416 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pt ! tracer field 417 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: ptc ! tracer content field 418 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(inout), OPTIONAL :: pua_ups, pva_ups ! upstream u*a fluxes 419 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out), OPTIONAL :: pua_ho, pva_ho ! high order u*a fluxes 273 420 ! 274 421 INTEGER :: ji, jj, jl ! dummy loop indices … … 303 450 DO jj = 1, jpjm1 304 451 DO ji = 1, fs_jpim1 305 IF( ABS( pu c(ji,jj,jl) ) > 0._wp .AND. ABS( pu(ji,jj) ) > 0._wp) THEN306 zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc (ji,jj,jl) / pu(ji,jj)307 zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * pu c(ji,jj,jl) / pu(ji,jj)452 IF( ABS( pu(ji,jj) ) > epsi10 ) THEN 453 zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc (ji,jj,jl) / pu(ji,jj) 454 zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * pua_ups(ji,jj,jl) / pu(ji,jj) 308 455 ELSE 309 456 zfu_ho (ji,jj,jl) = 0._wp … … 311 458 ENDIF 312 459 ! 313 IF( ABS( pv c(ji,jj,jl) ) > 0._wp .AND. ABS( pv(ji,jj) ) > 0._wp) THEN314 zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc (ji,jj,jl) / pv(ji,jj)315 zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pv c(ji,jj,jl) / pv(ji,jj)460 IF( ABS( pv(ji,jj) ) > epsi10 ) THEN 461 zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc (ji,jj,jl) / pv(ji,jj) 462 zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pva_ups(ji,jj,jl) / pv(ji,jj) 316 463 ELSE 317 464 zfv_ho (ji,jj,jl) = 0._wp … … 321 468 END DO 322 469 END DO 470 471 ! the new "volume" fluxes must also be "flux corrected" 472 ! thus we calculate the upstream solution and apply a limiter again 473 DO jl = 1, jpl 474 DO jj = 2, jpjm1 475 DO ji = fs_2, fs_jpim1 476 ztra = - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) + zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) 477 ! 478 zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1) 479 END DO 480 END DO 481 END DO 482 CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T', 1. ) 483 ! 484 IF ( np_limiter == 1 ) THEN 485 CALL nonosc_ice( 1._wp, pdt, pu, pv, ptc, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho ) 486 ELSEIF( np_limiter == 2 .OR. np_limiter == 3 ) THEN 487 CALL limiter_x( pdt, pu, ptc, zfu_ups, zfu_ho ) 488 CALL limiter_y( pdt, pv, ptc, zfv_ups, zfv_ho ) 489 ENDIF 490 ! 323 491 ENDIF 324 ! --ho 325 ! in case of advection of A: output u*a 326 ! ------------------------------------- 492 ! --ho --ups 493 ! in case of advection of A: output u*a and u*a 494 ! ----------------------------------------------- 327 495 IF( PRESENT( pua_ho ) ) THEN 328 496 DO jl = 1, jpl 329 497 DO jj = 1, jpjm1 330 498 DO ji = 1, fs_jpim1 331 pua_ho (ji,jj,jl) = zfu_ho(ji,jj,jl)332 p va_ho(ji,jj,jl) = zfv_ho(ji,jj,jl)333 499 pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 500 pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) 501 END DO 334 502 END DO 335 503 END DO … … 499 667 END DO 500 668 ! 501 IF ( kn_limiter == 1 ) THEN669 IF ( np_limiter == 1 ) THEN 502 670 CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 503 ELSEIF( kn_limiter == 2 .OR. kn_limiter == 3 ) THEN671 ELSEIF( np_limiter == 2 .OR. np_limiter == 3 ) THEN 504 672 CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 505 673 CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) … … 517 685 END DO 518 686 END DO 519 IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )687 IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 520 688 521 689 DO jl = 1, jpl !-- first guess of tracer from u-flux … … 538 706 END DO 539 707 END DO 540 IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )708 IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 541 709 542 710 ELSE !== even ice time step: adv_y then adv_x ==! … … 549 717 END DO 550 718 END DO 551 IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )719 IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 552 720 ! 553 721 DO jl = 1, jpl !-- first guess of tracer from v-flux … … 570 738 END DO 571 739 END DO 572 IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )740 IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 573 741 574 742 ENDIF 575 IF( kn_limiter == 1 ) CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )743 IF( np_limiter == 1 ) CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 576 744 577 745 ENDIF … … 609 777 ! 610 778 ! !-- ultimate interpolation of pt at u-point --! 611 CALL ultimate_x( kn_umx, pdt, pt, pu, zt_u, pfu_ho )779 CALL ultimate_x( pamsk, kn_umx, pdt, pt, pu, zt_u, pfu_ho ) 612 780 ! !-- limiter in x --! 613 IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )781 IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 614 782 ! !-- advective form update in zpt --! 615 783 DO jl = 1, jpl … … 619 787 & + pt (ji,jj,jl) * ( pu (ji,jj ) - pu (ji-1,jj ) ) * r1_e1e2t(ji,jj) & 620 788 & * pamsk & 621 & ) * pdt ) * tmask(ji,jj,1) 789 & ) * pdt ) * tmask(ji,jj,1) 622 790 END DO 623 791 END DO … … 627 795 ! !-- ultimate interpolation of pt at v-point --! 628 796 IF( ll_hoxy ) THEN 629 CALL ultimate_y( kn_umx, pdt, zpt, pv, zt_v, pfv_ho )797 CALL ultimate_y( pamsk, kn_umx, pdt, zpt, pv, zt_v, pfv_ho ) 630 798 ELSE 631 CALL ultimate_y( kn_umx, pdt, pt , pv, zt_v, pfv_ho )799 CALL ultimate_y( pamsk, kn_umx, pdt, pt , pv, zt_v, pfv_ho ) 632 800 ENDIF 633 801 ! !-- limiter in y --! 634 IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )802 IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 635 803 ! 636 804 ! … … 638 806 ! 639 807 ! !-- ultimate interpolation of pt at v-point --! 640 CALL ultimate_y( kn_umx, pdt, pt, pv, zt_v, pfv_ho )808 CALL ultimate_y( pamsk, kn_umx, pdt, pt, pv, zt_v, pfv_ho ) 641 809 ! !-- limiter in y --! 642 IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho )810 IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 643 811 ! !-- advective form update in zpt --! 644 812 DO jl = 1, jpl … … 656 824 ! !-- ultimate interpolation of pt at u-point --! 657 825 IF( ll_hoxy ) THEN 658 CALL ultimate_x( kn_umx, pdt, zpt, pu, zt_u, pfu_ho )826 CALL ultimate_x( pamsk, kn_umx, pdt, zpt, pu, zt_u, pfu_ho ) 659 827 ELSE 660 CALL ultimate_x( kn_umx, pdt, pt , pu, zt_u, pfu_ho )828 CALL ultimate_x( pamsk, kn_umx, pdt, pt , pu, zt_u, pfu_ho ) 661 829 ENDIF 662 830 ! !-- limiter in x --! 663 IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho )831 IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 664 832 ! 665 833 ENDIF 666 834 667 IF( kn_limiter == 1 ) CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )835 IF( np_limiter == 1 ) CALL nonosc_ice( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 668 836 ! 669 837 END SUBROUTINE macho 670 838 671 839 672 SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, pt_u, pfu_ho )840 SUBROUTINE ultimate_x( pamsk, kn_umx, pdt, pt, pu, pt_u, pfu_ho ) 673 841 !!--------------------------------------------------------------------- 674 842 !! *** ROUTINE ultimate_x *** … … 680 848 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 681 849 !!---------------------------------------------------------------------- 850 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) 682 851 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 683 852 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step … … 806 975 DO jj = 1, jpjm1 807 976 DO ji = 1, fs_jpim1 808 IF( pt_u(ji,jj,jl) < 0._wp ) THEN977 IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 809 978 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & 810 979 & - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) … … 826 995 827 996 828 SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pt_v, pfv_ho )997 SUBROUTINE ultimate_y( pamsk, kn_umx, pdt, pt, pv, pt_v, pfv_ho ) 829 998 !!--------------------------------------------------------------------- 830 999 !! *** ROUTINE ultimate_y *** … … 836 1005 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 837 1006 !!---------------------------------------------------------------------- 1007 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) 838 1008 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 839 1009 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step … … 959 1129 DO jj = 1, jpjm1 960 1130 DO ji = 1, fs_jpim1 961 IF( pt_v(ji,jj,jl) < 0._wp ) THEN1131 IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 962 1132 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) & 963 1133 & - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) … … 1023 1193 ! | | | | * 1024 1194 ! t_ups : i-1 i i+1 i+2 1025 IF( ll_prelim iter_zalesak) THEN1195 IF( ll_prelim ) THEN 1026 1196 1027 1197 DO jl = 1, jpl … … 1102 1272 ! 1103 1273 ! ! up & down beta terms 1104 IF( zpos > 0._wp ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt 1105 ELSE ; zbetup(ji,jj,jl) = 0._wp ! zbig 1274 ! clem: zbetup and zbetdo must be 0 for zpos>1.e-10 & zneg>1.e-10 (do not put 0 instead of 1.e-10 !!!) 1275 IF( zpos > epsi10 ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt 1276 ELSE ; zbetup(ji,jj,jl) = 0._wp ! zbig 1106 1277 ENDIF 1107 1278 ! 1108 IF( zneg > 0._wp) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt1109 ELSE ; zbetdo(ji,jj,jl) = 0._wp ! zbig1279 IF( zneg > epsi10 ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt 1280 ELSE ; zbetdo(ji,jj,jl) = 0._wp ! zbig 1110 1281 ENDIF 1111 1282 ! … … 1149 1320 END DO 1150 1321 1151 ! clem test1152 !! DO jj = 2, jpjm11153 !! DO ji = 2, fs_jpim1 ! vector opt.1154 !! zzt = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) &1155 !! & - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) &1156 !! & + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &1157 !! & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &1158 !! & ) * tmask(ji,jj,1)1159 !! IF( zzt < -epsi20 ) THEN1160 !! WRITE(numout,*) 'T<0 nonosc_ice',zzt1161 !! ENDIF1162 !! END DO1163 !! END DO1164 1165 1322 END DO 1166 1323 ! … … 1203 1360 Rjp = zslpx(ji+1,jj,jl) 1204 1361 1205 IF( kn_limiter == 3 ) THEN1362 IF( np_limiter == 3 ) THEN 1206 1363 1207 1364 IF( pu(ji,jj) > 0. ) THEN ; Rr = Rjm … … 1219 1376 pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter 1220 1377 1221 ELSEIF( kn_limiter == 2 ) THEN1378 ELSEIF( np_limiter == 2 ) THEN 1222 1379 IF( Rj /= 0. ) THEN 1223 1380 IF( pu(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj … … 1298 1455 Rjp = zslpy(ji,jj+1,jl) 1299 1456 1300 IF( kn_limiter == 3 ) THEN1457 IF( np_limiter == 3 ) THEN 1301 1458 1302 1459 IF( pv(ji,jj) > 0. ) THEN ; Rr = Rjm … … 1314 1471 pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter 1315 1472 1316 ELSEIF( kn_limiter == 2 ) THEN1473 ELSEIF( np_limiter == 2 ) THEN 1317 1474 1318 1475 IF( Rj /= 0. ) THEN … … 1358 1515 END SUBROUTINE limiter_y 1359 1516 1517 1518 SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 1519 !!------------------------------------------------------------------- 1520 !! *** ROUTINE Hbig *** 1521 !! 1522 !! ** Purpose : Thickness correction in case advection scheme creates 1523 !! abnormally tick ice or snow 1524 !! 1525 !! ** Method : 1- check whether ice thickness is larger than the surrounding 9-points 1526 !! (before advection) and reduce it by adapting ice concentration 1527 !! 2- check whether snow thickness is larger than the surrounding 9-points 1528 !! (before advection) and reduce it by sending the excess in the ocean 1529 !! 3- check whether snow load deplets the snow-ice interface below sea level$ 1530 !! and reduce it by sending the excess in the ocean 1531 !! 4- correct pond fraction to avoid a_ip > a_i 1532 !! 1533 !! ** input : Max thickness of the surrounding 9-points 1534 !!------------------------------------------------------------------- 1535 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 1536 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: phi_max, phs_max, phip_max ! max ice thick from surrounding 9-pts 1537 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip 1538 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s 1539 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i 1540 ! 1541 INTEGER :: ji, jj, jk, jl ! dummy loop indices 1542 REAL(wp) :: z1_dt, zhip, zhi, zhs, zvs_excess, zfra 1543 REAL(wp), DIMENSION(jpi,jpj) :: zswitch 1544 !!------------------------------------------------------------------- 1545 ! 1546 z1_dt = 1._wp / pdt 1547 ! 1548 DO jl = 1, jpl 1549 1550 DO jj = 1, jpj 1551 DO ji = 1, jpi 1552 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 1553 ! 1554 ! ! -- check h_ip -- ! 1555 ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 1556 IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 1557 zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 1558 IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 1559 pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl) 1560 ENDIF 1561 ENDIF 1562 ! 1563 ! ! -- check h_i -- ! 1564 ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i 1565 zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl) 1566 IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 1567 pa_i(ji,jj,jl) = pv_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) ) !-- bound h_i to hi_max (99 m) 1568 ENDIF 1569 ! 1570 ! ! -- check h_s -- ! 1571 ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean 1572 zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl) 1573 IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 1574 zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 ) 1575 ! 1576 wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * z1_dt 1577 hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 1578 ! 1579 pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 1580 pv_s(ji,jj,jl) = pa_i(ji,jj,jl) * phs_max(ji,jj,jl) 1581 ENDIF 1582 ! 1583 ! ! -- check snow load -- ! 1584 ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean 1585 ! this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin) 1586 ! this imposed mini can artificially make the snow very thick (if concentration decreases drastically) 1587 zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 1588 IF( zvs_excess > 0._wp ) THEN 1589 zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 ) 1590 wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt 1591 hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 1592 ! 1593 pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 1594 pv_s(ji,jj,jl) = pv_s(ji,jj,jl) - zvs_excess 1595 ENDIF 1596 1597 ENDIF 1598 END DO 1599 END DO 1600 END DO 1601 ! !-- correct pond fraction to avoid a_ip > a_i 1602 WHERE( pa_ip(:,:,:) > pa_i(:,:,:) ) pa_ip(:,:,:) = pa_i(:,:,:) 1603 ! 1604 ! 1605 END SUBROUTINE Hbig 1606 1360 1607 #else 1361 1608 !!---------------------------------------------------------------------- -
NEMO/branches/UKMO/NEMO_4.0_mirror/src/ICE/icedyn_rdgrft.F90
r10888 r11081 142 142 INTEGER, PARAMETER :: jp_itermax = 20 143 143 !!------------------------------------------------------------------- 144 ! clem: The redistribution of ice between categories can lead to small negative values (as for the remapping in ice_itd_rem)145 ! likely due to truncation error ( i.e. 1. - 1. /= 0 )146 ! I do not think it should be a concern since small areas and volumes are erased (in ice_var_zapsmall.F90)147 148 144 ! controls 149 145 IF( ln_timing ) CALL timing_start('icedyn_rdgrft') ! timing … … 156 152 ENDIF 157 153 158 CALL ice_var_zapsmall ! Zero out categories with very small areas159 160 154 !-------------------------------- 161 155 ! 0) Identify grid cells with ice 162 156 !-------------------------------- 157 at_i(:,:) = SUM( a_i, dim=3 ) 158 ! 163 159 npti = 0 ; nptidx(:) = 0 164 160 ipti = 0 ; iptidx(:) = 0 165 161 DO jj = 1, jpj 166 162 DO ji = 1, jpi 167 IF ( at_i(ji,jj) > 0._wp) THEN163 IF ( at_i(ji,jj) > epsi10 ) THEN 168 164 npti = npti + 1 169 165 nptidx( npti ) = (jj - 1) * jpi + ji … … 178 174 179 175 ! just needed here 180 CALL tab_2d_1d( npti, nptidx(1:npti), zdivu (1:npti), divu_i(:,:))181 CALL tab_2d_1d( npti, nptidx(1:npti), zdelt (1:npti), delta_i(:,:))176 CALL tab_2d_1d( npti, nptidx(1:npti), zdivu (1:npti) , divu_i ) 177 CALL tab_2d_1d( npti, nptidx(1:npti), zdelt (1:npti) , delta_i ) 182 178 ! needed here and in the iteration loop 183 CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:))184 CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:))185 CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti) , ato_i (:,:))179 CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i ) 180 CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i ) 181 CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti) , ato_i ) 186 182 187 183 DO ji = 1, npti … … 310 306 311 307 ! ! Ice thickness needed for rafting 312 WHERE( pa_i(1:npti,:) > 0._wp) ; zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)313 ELSEWHERE ; zhi(1:npti,:) = 0._wp308 WHERE( pa_i(1:npti,:) > epsi20 ) ; zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:) 309 ELSEWHERE ; zhi(1:npti,:) = 0._wp 314 310 END WHERE 315 311 … … 329 325 zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 ) 330 326 ! 331 WHERE( zasum(1:npti) > 0._wp) ; z1_asum(1:npti) = 1._wp / zasum(1:npti)332 ELSEWHERE ; z1_asum(1:npti) = 0._wp327 WHERE( zasum(1:npti) > epsi20 ) ; z1_asum(1:npti) = 1._wp / zasum(1:npti) 328 ELSEWHERE ; z1_asum(1:npti) = 0._wp 333 329 END WHERE 334 330 ! … … 455 451 ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate. 456 452 ! NOTE: 0 < aksum <= 1 457 WHERE( zaksum(1:npti) > 0._wp) ; closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)458 ELSEWHERE ; closing_gross(1:npti) = 0._wp453 WHERE( zaksum(1:npti) > epsi20 ) ; closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti) 454 ELSEWHERE ; closing_gross(1:npti) = 0._wp 459 455 END WHERE 460 456 … … 466 462 DO ji = 1, npti 467 463 zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice 468 IF( zfac > pa_i(ji,jl) ) THEN464 IF( zfac > pa_i(ji,jl) .AND. apartf(ji,jl) /= 0._wp ) THEN 469 465 closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_rdtice 470 466 ENDIF … … 510 506 REAL(wp), DIMENSION(jpij) :: zswitch, fvol ! new ridge volume going to jl2 511 507 REAL(wp), DIMENSION(jpij) :: z1_ai ! 1 / a 508 REAL(wp), DIMENSION(jpij) :: zvti ! sum(v_i) 512 509 ! 513 510 REAL(wp), DIMENSION(jpij,nlay_s) :: esrft ! snow energy of rafting ice … … 518 515 INTEGER , DIMENSION(jpij) :: itest_rdg, itest_rft ! test for conservation 519 516 !!------------------------------------------------------------------- 520 517 ! 518 zvti(1:npti) = SUM( v_i_2d(1:npti,:), dim=2 ) ! total ice volume 519 ! 521 520 ! 1) Change in open water area due to closing and opening 522 521 !-------------------------------------------------------- … … 535 534 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN ! only if ice is ridging 536 535 537 z1_ai(ji) = 1._wp / a_i_2d(ji,jl1) 538 536 IF( a_i_2d(ji,jl1) > epsi20 ) THEN ; z1_ai(ji) = 1._wp / a_i_2d(ji,jl1) 537 ELSE ; z1_ai(ji) = 0._wp 538 ENDIF 539 539 540 ! area of ridging / rafting ice (airdg1) and of new ridge (airdg2) 540 541 airdg1 = aridge(ji,jl1) * closing_gross(ji) * rdt_ice … … 549 550 550 551 ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges 551 vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg 552 IF ( zvti(ji) <= 10. ) THEN ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg ! v <= 10m then porosity = rn_porordg 553 ELSEIF( zvti(ji) >= 20. ) THEN ; vsw = 0._wp ! v >= 20m then porosity = 0 554 ELSE ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg * MAX( 0._wp, 2._wp - 0.1_wp * zvti(ji) ) ! v > 10m and v < 20m then porosity = linear transition to 0 555 ENDIF 552 556 ersw(ji) = -rhoi * vsw * rcp * sst_1d(ji) ! clem: if sst>0, then ersw <0 (is that possible?) 553 557 554 558 ! volume etc of ridging / rafting ice and new ridges (vi, vs, sm, oi, es, ei) 555 559 virdg1 = v_i_2d (ji,jl1) * afrdg 556 virdg2(ji) = v_i_2d (ji,jl1) * afrdg * ( 1. + rn_porordg )560 virdg2(ji) = v_i_2d (ji,jl1) * afrdg + vsw 557 561 vsrdg(ji) = v_s_2d (ji,jl1) * afrdg 558 562 sirdg1 = sv_i_2d(ji,jl1) * afrdg … … 726 730 END DO ! jl1 727 731 ! 732 ! roundoff errors 733 !---------------- 728 734 ! In case ridging/rafting lead to very small negative values (sometimes it happens) 729 WHERE( a_i_2d(1:npti,:) < 0._wp ) a_i_2d(1:npti,:) = 0._wp 730 WHERE( v_i_2d(1:npti,:) < 0._wp ) v_i_2d(1:npti,:) = 0._wp 735 CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d ) 731 736 ! 732 737 END SUBROUTINE rdgrft_shift … … 854 859 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 855 860 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd (:,:) ) 856 861 ! 857 862 ! !---------------------! 858 863 CASE( 2 ) !== from 1D to 2D ==! … … 945 950 CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' ) 946 951 ENDIF 947 ! ! allocate tke arrays 952 ! 953 IF( .NOT. ln_icethd ) THEN 954 rn_porordg = 0._wp 955 rn_fsnwrdg = 1._wp ; rn_fsnwrft = 1._wp 956 rn_fpndrdg = 1._wp ; rn_fpndrft = 1._wp 957 IF( lwp ) THEN 958 WRITE(numout,*) ' ==> only ice dynamics is activated, thus some parameters must be changed' 959 WRITE(numout,*) ' rn_porordg = ', rn_porordg 960 WRITE(numout,*) ' rn_fsnwrdg = ', rn_fsnwrdg 961 WRITE(numout,*) ' rn_fpndrdg = ', rn_fpndrdg 962 WRITE(numout,*) ' rn_fsnwrft = ', rn_fsnwrft 963 WRITE(numout,*) ' rn_fpndrft = ', rn_fpndrft 964 ENDIF 965 ENDIF 966 ! ! allocate arrays 948 967 IF( ice_dyn_rdgrft_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'ice_dyn_rdgrft_init: unable to allocate arrays' ) 949 968 ! -
NEMO/branches/UKMO/NEMO_4.0_mirror/src/ICE/icedyn_rhg.F90
r10888 r11081 58 58 !!-------------------------------------------------------------------- 59 59 INTEGER, INTENT(in) :: kt ! ice time step 60 ! 61 INTEGER :: jl ! dummy loop indices 60 62 !!-------------------------------------------------------------------- 61 63 ! controls … … 68 70 WRITE(numout,*)'~~~~~~~~~~~' 69 71 ENDIF 70 71 ! -------- 72 ! Rheology 73 ! -------- 72 ! 73 IF( ln_landfast_home ) THEN !-- Landfast ice parameterization 74 tau_icebfr(:,:) = 0._wp 75 DO jl = 1, jpl 76 WHERE( h_i(:,:,jl) > ht_n(:,:) * rn_depfra ) tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 77 END DO 78 ENDIF 79 ! 80 !--------------! 81 !== Rheology ==! 82 !--------------! 74 83 SELECT CASE( nice_rhg ) 75 84 ! !------------------------! -
NEMO/branches/UKMO/NEMO_4.0_mirror/src/ICE/icedyn_rhg_evp.F90
r10888 r11081 112 112 REAL(wp), DIMENSION(:,:), INTENT( out) :: pshear_i , pdivu_i , pdelta_i ! 113 113 !! 114 LOGICAL, PARAMETER :: ll_bdy_substep = . FALSE. ! temporary option to call bdy at each sub-time step (T)114 LOGICAL, PARAMETER :: ll_bdy_substep = .TRUE. ! temporary option to call bdy at each sub-time step (T) 115 115 ! or only at the main time step (F) 116 116 INTEGER :: ji, jj ! dummy loop indices … … 160 160 161 161 REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter 162 REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity equals ocean velocity 162 REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity becomes very small 163 REAL(wp), PARAMETER :: zamin = 0.001_wp ! ice concentration below which ice velocity becomes very small 163 164 !! --- diags 164 165 REAL(wp), DIMENSION(jpi,jpj) :: zswi … … 319 320 320 321 ! switches 321 zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin 322 zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin 322 IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN ; zswitchU(ji,jj) = 0._wp 323 ELSE ; zswitchU(ji,jj) = 1._wp ; ENDIF 324 IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN ; zswitchV(ji,jj) = 0._wp 325 ELSE ; zswitchV(ji,jj) = 1._wp ; ENDIF 323 326 324 327 END DO … … 519 522 & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 520 523 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 521 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin524 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 522 525 & ) * zmaskV(ji,jj) 523 526 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) … … 526 529 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 527 530 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 528 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin531 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 529 532 & ) * zmaskV(ji,jj) 530 533 ENDIF … … 567 570 & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 568 571 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 569 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin572 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 570 573 & ) * zmaskU(ji,jj) 571 574 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) … … 574 577 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 575 578 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 576 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin579 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 577 580 & ) * zmaskU(ji,jj) 578 581 ENDIF … … 617 620 & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 618 621 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 619 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin622 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 620 623 & ) * zmaskU(ji,jj) 621 624 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) … … 624 627 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 625 628 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 626 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin629 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 627 630 & ) * zmaskU(ji,jj) 628 631 ENDIF … … 665 668 & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 666 669 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 667 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin670 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 668 671 & ) * zmaskV(ji,jj) 669 672 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) … … 672 675 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 673 676 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 674 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin677 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 675 678 & ) * zmaskV(ji,jj) 676 679 ENDIF -
NEMO/branches/UKMO/NEMO_4.0_mirror/src/ICE/iceitd.F90
r10888 r11081 21 21 USE ice1D ! sea-ice: thermodynamic variables 22 22 USE ice ! sea-ice: variables 23 USE icevar ! sea-ice: operations 23 24 USE icectl ! sea-ice: conservation tests 24 25 USE icetab ! sea-ice: convert 1D<=>2D … … 91 92 ! 1) Identify grid cells with ice 92 93 !----------------------------------------------------------------------------------------------- 94 at_i(:,:) = SUM( a_i, dim=3 ) 95 ! 93 96 npti = 0 ; nptidx(:) = 0 94 97 DO jj = 1, jpj … … 249 252 ! --- g(h) for each thickness category --- ! 250 253 CALL itd_glinear( zhbnew(1:npti,jl-1), zhbnew(1:npti,jl), h_i_1d(1:npti) , a_i_1d(1:npti) , & ! in 251 & g0 (1:npti,jl ), g1 (1:npti,jl), hL (1:npti,jl), hR(1:npti,jl) ) ! out254 & g0 (1:npti,jl ), g1 (1:npti,jl), hL (1:npti,jl), hR (1:npti,jl) ) ! out 252 255 ! 253 256 END DO … … 389 392 REAL(wp), DIMENSION(:,:), INTENT(in) :: pdvice ! ice volume transferred across boundary 390 393 ! 391 INTEGER :: ji, j j, jl, jk! dummy loop indices392 INTEGER :: ii, ij, jl2, jl1! local integers394 INTEGER :: ji, jl, jk ! dummy loop indices 395 INTEGER :: jl2, jl1 ! local integers 393 396 REAL(wp) :: ztrans ! ice/snow transferred 394 REAL(wp), DIMENSION(jpij) :: zworka, zworkv ! workspace 395 REAL(wp), DIMENSION(jpij,jpl) :: zaTsfn ! - - 397 REAL(wp), DIMENSION(jpij) :: zworka, zworkv ! workspace 398 REAL(wp), DIMENSION(jpij,jpl) :: zaTsfn ! - - 399 REAL(wp), DIMENSION(jpij,nlay_i,jpl) :: ze_i_2d 400 REAL(wp), DIMENSION(jpij,nlay_s,jpl) :: ze_s_2d 396 401 !!------------------------------------------------------------------ 397 402 … … 405 410 CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip ) 406 411 CALL tab_3d_2d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su ) 412 DO jl = 1, jpl 413 DO jk = 1, nlay_s 414 CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) ) 415 END DO 416 DO jk = 1, nlay_i 417 CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) ) 418 END DO 419 END DO 420 ! to correct roundoff errors on a_i 421 CALL tab_2d_1d( npti, nptidx(1:npti), rn_amax_1d(1:npti), rn_amax_2d ) 407 422 408 423 !---------------------------------------------------------------------------------------------- … … 435 450 ELSE ; zworka(ji) = 0._wp 436 451 ENDIF 437 !438 ! clem: The transfer between one category to another can lead to very small negative values (-1.e-20)439 ! because of truncation error ( i.e. 1. - 1. /= 0 )440 ! I do not think it should be a concern since small areas and volumes are erased (in ice_var_zapsmall.F90)441 452 ! 442 453 a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl) ! Ice areas … … 476 487 ! 477 488 DO jk = 1, nlay_s !--- Snow heat content 478 !479 489 DO ji = 1, npti 480 ii = MOD( nptidx(ji) - 1, jpi ) + 1481 ij = ( nptidx(ji) - 1 ) / jpi + 1482 490 ! 483 491 jl1 = kdonor(ji,jl) … … 487 495 ELSE ; jl2 = jl 488 496 ENDIF 489 ! 490 ztrans = e_s(ii,ij,jk,jl1) * zworkv(ji) 491 e_s(ii,ij,jk,jl1) = e_s(ii,ij,jk,jl1) - ztrans 492 e_s(ii,ij,jk,jl2) = e_s(ii,ij,jk,jl2) + ztrans 497 ztrans = ze_s_2d(ji,jk,jl1) * zworkv(ji) 498 ze_s_2d(ji,jk,jl1) = ze_s_2d(ji,jk,jl1) - ztrans 499 ze_s_2d(ji,jk,jl2) = ze_s_2d(ji,jk,jl2) + ztrans 493 500 ENDIF 494 501 END DO … … 497 504 DO jk = 1, nlay_i !--- Ice heat content 498 505 DO ji = 1, npti 499 ii = MOD( nptidx(ji) - 1, jpi ) + 1500 ij = ( nptidx(ji) - 1 ) / jpi + 1501 506 ! 502 507 jl1 = kdonor(ji,jl) … … 506 511 ELSE ; jl2 = jl 507 512 ENDIF 508 ! 509 ztrans = e_i(ii,ij,jk,jl1) * zworkv(ji) 510 e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - ztrans 511 e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + ztrans 513 ztrans = ze_i_2d(ji,jk,jl1) * zworkv(ji) 514 ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) - ztrans 515 ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + ztrans 512 516 ENDIF 513 517 END DO … … 515 519 ! 516 520 END DO ! boundaries, 1 to jpl-1 521 522 !------------------- 523 ! 3) roundoff errors 524 !------------------- 525 ! clem: The transfer between one category to another can lead to very small negative values (-1.e-20) 526 ! because of truncation error ( i.e. 1. - 1. /= 0 ) 527 CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d ) 528 529 ! at_i must be <= rn_amax 530 zworka(1:npti) = SUM( a_i_2d(1:npti,:), dim=2 ) 531 DO jl = 1, jpl 532 WHERE( zworka(1:npti) > rn_amax_1d(1:npti) ) & 533 & a_i_2d(1:npti,jl) = a_i_2d(1:npti,jl) * rn_amax_1d(1:npti) / zworka(1:npti) 534 END DO 517 535 518 536 !------------------------------------------------------------------------------- 519 ! 3) Update ice thickness and temperature537 ! 4) Update ice thickness and temperature 520 538 !------------------------------------------------------------------------------- 521 539 WHERE( a_i_2d(1:npti,:) >= epsi20 ) … … 536 554 CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip ) 537 555 CALL tab_2d_3d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su ) 556 DO jl = 1, jpl 557 DO jk = 1, nlay_s 558 CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) ) 559 END DO 560 DO jk = 1, nlay_i 561 CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) ) 562 END DO 563 END DO 538 564 ! 539 565 END SUBROUTINE itd_shiftice … … 558 584 ! 559 585 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_itd_reb: rebining ice thickness distribution' 586 ! 587 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 560 588 ! 561 589 jdonor(:,:) = 0 … … 635 663 END DO 636 664 ! 665 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 666 ! 637 667 END SUBROUTINE ice_itd_reb 638 668 -
NEMO/branches/UKMO/NEMO_4.0_mirror/src/ICE/icestp.F90
r10888 r11081 189 189 IF( ln_icethd ) CALL ice_thd( kt ) ! -- Ice thermodynamics 190 190 ! 191 IF( ln_icethd )CALL ice_cor( kt , 2 ) ! -- Corrections191 CALL ice_cor( kt , 2 ) ! -- Corrections 192 192 ! 193 193 CALL ice_var_glo2eqv ! necessary calls (at least for coupling) … … 323 323 WRITE(numout,*) ' maximum ice concentration for SH = ', rn_amax_s 324 324 ENDIF 325 ! !--- change max ice concentration for roundoff errors 326 rn_amax_n = MIN( rn_amax_n, 1._wp - epsi10 ) 327 rn_amax_s = MIN( rn_amax_s, 1._wp - epsi10 ) 325 328 ! !--- check consistency 326 329 IF ( jpl > 1 .AND. ln_virtual_itd ) THEN … … 431 434 t_si (:,:,:) = rt0 ! temp at the ice-snow interface 432 435 433 tau_icebfr(:,:) = 0._wp ! landfast ice param only (clem: important to keep the init here) 434 cnd_ice (:,:,:) = 0._wp ! initialisation: effective conductivity at the top of ice/snow (ln_cndflx=T) 435 qtr_ice_bot(:,:,:) = 0._wp ! initialization: part of solar radiation transmitted through the ice needed at least for outputs 436 tau_icebfr (:,:) = 0._wp ! landfast ice param only (clem: important to keep the init here) 437 cnd_ice (:,:,:) = 0._wp ! initialisation: effective conductivity at the top of ice/snow (ln_cndflx=T) 438 qcn_ice (:,:,:) = 0._wp ! initialisation: conductive flux (ln_cndflx=T & ln_cndemule=T) 439 qtr_ice_bot(:,:,:) = 0._wp ! initialization: part of solar radiation transmitted through the ice needed at least for outputs 440 qsb_ice_bot(:,:) = 0._wp ! (needed if ln_icethd=F) 436 441 ! 437 442 ! for control checks (ln_icediachk) -
NEMO/branches/UKMO/NEMO_4.0_mirror/src/ICE/icethd.F90
r10888 r11081 102 102 ENDIF 103 103 104 CALL ice_var_glo2eqv105 106 104 !---------------------------------------------! 107 105 ! computation of friction velocity at T points … … 162 160 qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rdt_ice ) - zqfr ) 163 161 164 ! If there is ice and leads are warming, then transfer energy from the lead budget and use it for bottom melting 165 IF( zqld > 0._wp ) THEN 162 ! If there is ice and leads are warming => transfer energy from the lead budget and use it for bottom melting 163 ! If the grid cell is fully covered by ice (no leads) => transfer energy from the lead budget to the ice bottom budget 164 IF( ( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) .OR. at_i(ji,jj) >= (1._wp - epsi10) ) THEN 166 165 fhld (ji,jj) = rswitch * zqld * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90 167 166 qlead(ji,jj) = 0._wp … … 178 177 ! In case we bypass open-water ice formation 179 178 IF( .NOT. ln_icedO ) qlead(:,:) = 0._wp 180 ! In case we bypass growing/melting from top and bottom : we suppose ice is impermeable => ocean is isolated from atmosphere179 ! In case we bypass growing/melting from top and bottom 181 180 IF( .NOT. ln_icedH ) THEN 182 qt_atm_oi (:,:) = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:)183 181 qsb_ice_bot(:,:) = 0._wp 184 182 fhld (:,:) = 0._wp … … 221 219 dh_i_sub (1:npti) = 0._wp ; dh_i_bog(1:npti) = 0._wp 222 220 dh_snowice(1:npti) = 0._wp ; dh_s_mlt(1:npti) = 0._wp 223 ! 224 IF( ln_icedH ) THEN ! --- growing/melting --- ! 225 CALL ice_thd_zdf ! Ice/Snow Temperature profile 226 CALL ice_thd_dh ! Ice/Snow thickness 227 CALL ice_thd_pnd ! Melt ponds formation 228 CALL ice_thd_ent( e_i_1d(1:npti,:) ) ! Ice enthalpy remapping 221 ! 222 CALL ice_thd_zdf ! --- Ice-Snow temperature --- ! 223 ! 224 IF( ln_icedH ) THEN ! --- Growing/Melting --- ! 225 CALL ice_thd_dh ! Ice-Snow thickness 226 CALL ice_thd_pnd ! Melt ponds formation 227 CALL ice_thd_ent( e_i_1d(1:npti,:) ) ! Ice enthalpy remapping 229 228 ENDIF 230 !231 229 CALL ice_thd_sal( ln_icedS ) ! --- Ice salinity --- ! 232 230 ! 233 CALL ice_thd_temp ! --- temperature update --- !231 CALL ice_thd_temp ! --- Temperature update --- ! 234 232 ! 235 233 IF( ln_icedH .AND. ln_virtual_itd ) & 236 & CALL ice_thd_mono ! --- extra lateral melting if virtual_itd --- !237 ! 238 IF( ln_icedA ) CALL ice_thd_da ! --- lateral melting --- !234 & CALL ice_thd_mono ! --- Extra lateral melting if virtual_itd --- ! 235 ! 236 IF( ln_icedA ) CALL ice_thd_da ! --- Lateral melting --- ! 239 237 ! 240 238 CALL ice_thd_1d2d( jl, 2 ) ! --- Change units of e_i, e_s from J/m3 to J/m2 --- ! 241 239 ! ! --- & Move to 2D arrays --- ! 242 !243 240 ENDIF 244 241 ! 245 242 END DO 246 243 ! 247 244 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 248 !249 CALL ice_var_zapsmall ! --- remove very small ice concentration (<1e-10) --- !250 ! ! & make sure at_i=SUM(a_i) & ato_i=1 where at_i=0251 245 ! 252 IF( jpl > 1 ) CALL ice_itd_rem( kt )! --- Transport ice between thickness categories --- !253 ! 254 IF( ln_icedO ) CALL ice_thd_do ! --- frazil ice growingin leads --- !246 IF( jpl > 1 ) CALL ice_itd_rem( kt ) ! --- Transport ice between thickness categories --- ! 247 ! 248 IF( ln_icedO ) CALL ice_thd_do ! --- Frazil ice growth in leads --- ! 255 249 ! 256 250 ! controls -
NEMO/branches/UKMO/NEMO_4.0_mirror/src/ICE/icethd_do.F90
r10888 r11081 114 114 IF( ln_icediachk ) CALL ice_cons_hsm( 0, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft ) 115 115 116 CALL ice_var_agg(1) 117 CALL ice_var_glo2eqv 118 116 at_i(:,:) = SUM( a_i, dim=3 ) 119 117 !------------------------------------------------------------------------------! 120 118 ! 1) Collection thickness of ice formed in leads and polynyas … … 319 317 320 318 ! --- lateral ice growth --- ! 321 ! If lateral ice growth gives an ice concentration gt 1, then319 ! If lateral ice growth gives an ice concentration > amax, then 322 320 ! we keep the excessive volume in memory and attribute it later to bottom accretion 323 321 DO ji = 1, npti 324 IF ( za_newice(ji) > ( rn_amax_1d(ji) - at_i_1d(ji) ) ) THEN325 zda_res(ji) = za_newice(ji) - (rn_amax_1d(ji) - at_i_1d(ji) )322 IF ( za_newice(ji) > MAX( 0._wp, rn_amax_1d(ji) - at_i_1d(ji) ) ) THEN ! max is for roundoff error 323 zda_res(ji) = za_newice(ji) - MAX( 0._wp, rn_amax_1d(ji) - at_i_1d(ji) ) 326 324 zdv_res(ji) = zda_res (ji) * zh_newice(ji) 327 za_newice(ji) = za_newice(ji) - zda_res (ji)328 zv_newice(ji) = zv_newice(ji) - zdv_res (ji)325 za_newice(ji) = MAX( 0._wp, za_newice(ji) - zda_res (ji) ) 326 zv_newice(ji) = MAX( 0._wp, zv_newice(ji) - zdv_res (ji) ) 329 327 ELSE 330 328 zda_res(ji) = 0._wp -
NEMO/branches/UKMO/NEMO_4.0_mirror/src/ICE/icethd_zdf_bl99.F90
r10888 r11081 206 206 ! 207 207 l_T_converged(:) = .FALSE. 208 ! !============================!209 208 ! Convergence calculated until all sub-domain grid points have converged 210 209 ! Calculations keep going for all grid points until sub-domain convergence (vectorisation optimisation) 211 210 ! but values are not taken into account (results independant of MPI partitioning) 212 211 ! 212 ! !============================! 213 213 DO WHILE ( ( .NOT. ALL (l_T_converged(1:npti)) ) .AND. iconv < iconv_max ) ! Iterative procedure begins ! 214 ! !============================!214 ! !============================! 215 215 iconv = iconv + 1 216 216 ! … … 742 742 zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) ) 743 743 ! t_i 744 DO jk = 0, nlay_i744 DO jk = 1, nlay_i 745 745 ztmelts = -rTmlt * sz_i_1d(ji,jk) + rt0 746 746 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) … … 856 856 t_i_1d (1:npti,:) = ztiold (1:npti,:) 857 857 qcn_ice_1d(1:npti) = qcn_ice_top_1d(1:npti) 858 859 !!clem 860 ! remettre t_su_1d, qns_ice_1d et dqns_ice_1d comme avant puisqu'on devrait faire comme si on avant conduction = input 861 !clem 858 862 ENDIF 859 863 ! -
NEMO/branches/UKMO/NEMO_4.0_mirror/src/ICE/icevar.F90
r10888 r11081 44 44 !! ice_var_salprof1d : salinity profile in the ice 1D 45 45 !! ice_var_zapsmall : remove very small area and volume 46 !! ice_var_zapneg : remove negative ice fields (to debug the advection scheme UM3-5) 46 !! ice_var_zapneg : remove negative ice fields 47 !! ice_var_roundoff : remove negative values arising from roundoff erros 47 48 !! ice_var_itd : convert 1-cat to jpl-cat 48 49 !! ice_var_itd2 : convert N-cat to jpl-cat … … 71 72 PUBLIC ice_var_zapsmall 72 73 PUBLIC ice_var_zapneg 74 PUBLIC ice_var_roundoff 73 75 PUBLIC ice_var_itd 74 76 PUBLIC ice_var_itd2 … … 229 231 IF ( v_i(ji,jj,jl) > epsi20 ) THEN !--- icy area 230 232 ! 231 ze_i = e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i 233 ze_i = e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i ! Energy of melting e(S,T) [J.m-3] 232 234 ztmelts = - sz_i(ji,jj,jk,jl) * rTmlt ! Ice layer melt temperature [C] 233 235 ! Conversion q(S,T) -> T (second order equation) … … 236 238 t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_rcpi , ztmelts ) ) + rt0 ! [K] with bounds: -100 < t_i < ztmelts 237 239 ! 238 ELSE !--- no ice240 ELSE !--- no ice 239 241 t_i(ji,jj,jk,jl) = rt0 240 242 ENDIF … … 537 539 538 540 539 SUBROUTINE ice_var_zapneg( p ato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )541 SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 540 542 !!------------------------------------------------------------------- 541 543 !! *** ROUTINE ice_var_zapneg *** … … 543 545 !! ** Purpose : Remove negative sea ice fields and correct fluxes 544 546 !!------------------------------------------------------------------- 545 INTEGER :: ji, jj, jl, jk ! dummy loop indices 546 ! 547 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 547 548 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pato_i ! open water area 548 549 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i ! ice volume … … 555 556 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content 556 557 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content 557 !!------------------------------------------------------------------- 558 ! 558 ! 559 INTEGER :: ji, jj, jl, jk ! dummy loop indices 560 REAL(wp) :: z1_dt 561 !!------------------------------------------------------------------- 562 ! 563 z1_dt = 1._wp / pdt 559 564 ! 560 565 DO jl = 1, jpl !== loop over the categories ==! 561 566 ! 567 ! make sure a_i=0 where v_i<=0 568 WHERE( pv_i(:,:,:) <= 0._wp ) pa_i(:,:,:) = 0._wp 569 562 570 !---------------------------------------- 563 571 ! zap ice energy and send it to the ocean … … 566 574 DO jj = 1 , jpj 567 575 DO ji = 1 , jpi 568 IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN569 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * r1_rdtice! W.m-2 >0576 IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 577 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0 570 578 pe_i(ji,jj,jk,jl) = 0._wp 571 579 ENDIF … … 577 585 DO jj = 1 , jpj 578 586 DO ji = 1 , jpi 579 IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN580 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * r1_rdtice! W.m-2 <0587 IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 588 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0 581 589 pe_s(ji,jj,jk,jl) = 0._wp 582 590 ENDIF … … 590 598 DO jj = 1 , jpj 591 599 DO ji = 1 , jpi 592 IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN593 wfx_res(ji,jj) = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * r1_rdtice600 IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 601 wfx_res(ji,jj) = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt 594 602 pv_i (ji,jj,jl) = 0._wp 595 603 ENDIF 596 IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN597 wfx_res(ji,jj) = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * r1_rdtice604 IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 605 wfx_res(ji,jj) = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt 598 606 pv_s (ji,jj,jl) = 0._wp 599 607 ENDIF 600 IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN601 sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * r1_rdtice608 IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 609 sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt 602 610 psv_i (ji,jj,jl) = 0._wp 603 611 ENDIF … … 616 624 END SUBROUTINE ice_var_zapneg 617 625 626 627 SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pe_s, pe_i ) 628 !!------------------------------------------------------------------- 629 !! *** ROUTINE ice_var_roundoff *** 630 !! 631 !! ** Purpose : Remove negative sea ice values arising from roundoff errors 632 !!------------------------------------------------------------------- 633 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pa_i ! ice concentration 634 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_i ! ice volume 635 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_s ! snw volume 636 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: psv_i ! salt content 637 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: poa_i ! age content 638 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pa_ip ! melt pond fraction 639 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_ip ! melt pond volume 640 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_s ! snw heat content 641 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_i ! ice heat content 642 !!------------------------------------------------------------------- 643 ! 644 WHERE( pa_i (1:npti,:) < 0._wp .AND. pa_i (1:npti,:) > -epsi10 ) pa_i (1:npti,:) = 0._wp ! a_i must be >= 0 645 WHERE( pv_i (1:npti,:) < 0._wp .AND. pv_i (1:npti,:) > -epsi10 ) pv_i (1:npti,:) = 0._wp ! v_i must be >= 0 646 WHERE( pv_s (1:npti,:) < 0._wp .AND. pv_s (1:npti,:) > -epsi10 ) pv_s (1:npti,:) = 0._wp ! v_s must be >= 0 647 WHERE( psv_i(1:npti,:) < 0._wp .AND. psv_i(1:npti,:) > -epsi10 ) psv_i(1:npti,:) = 0._wp ! sv_i must be >= 0 648 WHERE( poa_i(1:npti,:) < 0._wp .AND. poa_i(1:npti,:) > -epsi10 ) poa_i(1:npti,:) = 0._wp ! oa_i must be >= 0 649 WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 ) pe_i (1:npti,:,:) = 0._wp ! e_i must be >= 0 650 WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 ) pe_s (1:npti,:,:) = 0._wp ! e_s must be >= 0 651 IF ( ln_pnd_H12 ) THEN 652 WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 ) pa_ip(1:npti,:) = 0._wp ! a_ip must be >= 0 653 WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 ) pv_ip(1:npti,:) = 0._wp ! v_ip must be >= 0 654 ENDIF 655 ! 656 END SUBROUTINE ice_var_roundoff 657 618 658 619 659 SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i ) … … 650 690 INTEGER :: idim, i_fill, jl0 651 691 REAL(wp) :: zarg, zV, zconv, zdh, zdv 652 REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables653 REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables692 REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables 693 REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables 654 694 INTEGER , DIMENSION(4) :: itest 655 695 !!------------------------------------------------------------------- … … 780 820 !! 781 821 !! 2) Expand the filling to the cat jlmin-1 and jlmax+1 782 !! by removing 25% ice area from jlmin and jlmax (resp.)822 !! by removing 25% ice area from jlmin and jlmax (resp.) 783 823 !! 784 824 !! 3) Expand the filling to the empty cat between jlmin and jlmax -
NEMO/branches/UKMO/NEMO_4.0_mirror/src/ICE/icewri.F90
r10888 r11081 145 145 146 146 ! record presence of fast ice 147 WHERE( z2d(:,:) < 5.e-04_wp .AND. zmsk 00(:,:) == 1._wp ) ; zfast(:,:) = 1._wp147 WHERE( z2d(:,:) < 5.e-04_wp .AND. zmsk15(:,:) == 1._wp ) ; zfast(:,:) = 1._wp 148 148 ELSEWHERE ; zfast(:,:) = 0._wp 149 149 END WHERE
Note: See TracChangeset
for help on using the changeset viewer.