- Timestamp:
- 2019-06-06T16:11:54+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
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 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
