Changeset 10911 for NEMO/trunk/src/ICE/icedyn_adv_umx.F90
- Timestamp:
- 2019-04-29T18:11:12+02:00 (5 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/ICE/icedyn_adv_umx.F90 (modified) (35 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/ICE/icedyn_adv_umx.F90
r10786 r10911 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 … … 39 39 INTEGER :: kn_limiter = 1 40 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 41 ! if there is an outward velocity in a grid cell where there is no ice initially (typically at the ice edge), 42 ! interpolated T at u/v points can be non-zero while it should 43 ! (because of the high order of the advection scheme). Thus set it to 0 in this case 44 LOGICAL :: ll_icedge = .TRUE. 45 46 ! if T interpolated at u/v points is negative or v_i < 1.e-6 47 ! then interpolate T at u/v points using the upstream scheme 44 48 LOGICAL :: ll_neg = .TRUE. 45 49 … … 51 55 52 56 ! 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 57 LOGICAL :: ll_prelimiter_zalesak = .FALSE. ! from: Zalesak(1979) eq. 14 => better for 1D. Not well defined in 2D 55 58 59 ! advection for S and T: dVS/dt = -div( uA * uHS / u ) => ll_advS = F 60 ! or dVS/dt = -div( uV * uS / u ) => ll_advS = T 61 LOGICAL :: ll_advS = .FALSE. 56 62 57 63 !! * Substitutions … … 64 70 CONTAINS 65 71 66 SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, &72 SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip, & 67 73 & pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 68 74 !!---------------------------------------------------------------------- … … 79 85 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pu_ice ! ice i-velocity 80 86 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pv_ice ! ice j-velocity 87 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_i ! ice thickness 88 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_s ! snw thickness 89 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_ip ! ice pond thickness 81 90 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pato_i ! open water area 82 91 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i ! ice volume … … 94 103 REAL(wp) :: zamsk ! 1 if advection of concentration, 0 if advection of other tracers 95 104 REAL(wp) :: zdt 96 REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! send zcflnow and receive zcflprv105 REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! for global communication 97 106 REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx, zcu_box, zcv_box 98 107 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 108 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zua_ho, zva_ho, zua_ups, zva_ups 109 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zuv_ho, zvv_ho, zuv_ups, zvv_ups 110 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_ai, z1_aip, z1_vi, z1_vs 101 111 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhvar 112 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhi_max, zhs_max, zhip_max 102 113 !!---------------------------------------------------------------------- 103 114 ! 104 115 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme' 105 116 ! 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... 117 ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- ! 118 DO jl = 1, jpl 119 DO jj = 2, jpjm1 120 DO ji = fs_2, fs_jpim1 121 zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj ,jl), ph_ip(ji ,jj+1,jl), & 122 & ph_ip(ji-1,jj ,jl), ph_ip(ji ,jj-1,jl), & 123 & ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 124 & ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 125 zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj ,jl), ph_i (ji ,jj+1,jl), & 126 & ph_i (ji-1,jj ,jl), ph_i (ji ,jj-1,jl), & 127 & ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 128 & ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 129 zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj ,jl), ph_s (ji ,jj+1,jl), & 130 & ph_s (ji-1,jj ,jl), ph_s (ji ,jj-1,jl), & 131 & ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 132 & ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 133 END DO 134 END DO 135 END DO 136 CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1., zhip_max(:,:,:), 'T', 1. ) 137 ! 138 ! 139 ! --- If ice drift is too fast, use subtime steps for advection (CFL test for stability) --- ! 140 ! Note: the advection split is applied at the next time-step in order to avoid blocking global comm. 141 ! this should not affect too much the stability 109 142 zcflnow(1) = MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) 110 143 zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) … … 116 149 ELSE ; icycle = 1 117 150 ENDIF 118 119 151 zdt = rdt_ice / REAL(icycle) 120 152 … … 154 186 END WHERE 155 187 ! 156 ! set u* a=u for advection of A only188 ! set u*A=u for advection of A only 157 189 DO jl = 1, jpl 158 190 zua_ho(:,:,jl) = zudy(:,:) 159 191 zva_ho(:,:,jl) = zvdx(:,:) 160 192 END DO 161 193 194 !== Ice area ==! 162 195 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 196 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 197 & pa_i, pa_i, zua_ups, zva_ups, zua_ho , zva_ho ) 164 198 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 199 200 IF( .NOT. ll_advS ) THEN !-- advection form: uA * uHS / u 201 !== Ice volume ==! 202 zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:) 203 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 204 & zhvar, pv_i, zua_ups, zva_ups ) 205 ! 206 !== Snw volume ==! 207 zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:) 208 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 209 & zhvar, pv_s, zua_ups, zva_ups ) 210 ! 211 !== Salt content ==! 212 zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:) 213 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 214 & zhvar, psv_i, zua_ups, zva_ups ) 215 ! 216 !== Ice heat content ==! 217 DO jk = 1, nlay_i 218 zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:) 219 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho, zva_ho, zcu_box, zcv_box, & 220 & zhvar, pe_i(:,:,jk,:), zua_ups, zva_ups ) 221 END DO 222 ! 223 !== Snw heat content ==! 224 DO jk = 1, nlay_s 225 zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:) 226 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho, zva_ho, zcu_box, zcv_box, & 227 & zhvar, pe_s(:,:,jk,:), zua_ups, zva_ups ) 228 END DO 229 ! 230 ELSE !-- advection form: uV * uS / u 231 ! inverse of Vi 232 WHERE( pv_i(:,:,:) >= epsi20 ) ; z1_vi(:,:,:) = 1._wp / pv_i(:,:,:) 233 ELSEWHERE ; z1_vi(:,:,:) = 0. 234 END WHERE 235 ! inverse of Vs 236 WHERE( pv_s(:,:,:) >= epsi20 ) ; z1_vs(:,:,:) = 1._wp / pv_s(:,:,:) 237 ELSEWHERE ; z1_vs(:,:,:) = 0. 238 END WHERE 239 ! 240 ! It is important to first calculate the ice fields and then the snow fields (because we use the same arrays) 241 ! 242 !== Ice volume ==! 243 zuv_ups = zua_ups 244 zvv_ups = zva_ups 245 zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:) 246 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 247 & zhvar, pv_i, zuv_ups, zvv_ups, zuv_ho , zvv_ho ) 248 ! 249 !== Salt content ==! 250 zhvar(:,:,:) = psv_i(:,:,:) * z1_vi(:,:,:) 251 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zuv_ho , zvv_ho , zcu_box, zcv_box, & 252 & zhvar, psv_i, zuv_ups, zvv_ups ) 253 ! 254 !== Ice heat content ==! 255 DO jk = 1, nlay_i 256 zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_vi(:,:,:) 257 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zuv_ho, zvv_ho, zcu_box, zcv_box, & 258 & zhvar, pe_i(:,:,jk,:), zuv_ups, zvv_ups ) 259 END DO 260 ! 261 !== Snow volume ==! 262 zuv_ups = zua_ups 263 zvv_ups = zva_ups 264 zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:) 265 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 266 & zhvar, pv_s, zuv_ups, zvv_ups, zuv_ho , zvv_ho ) 267 ! 268 !== Snw heat content ==! 269 DO jk = 1, nlay_s 270 zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_vs(:,:,:) 271 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zuv_ho, zvv_ho, zcu_box, zcv_box, & 272 & zhvar, pe_s(:,:,jk,:), zuv_ups, zvv_ups ) 273 END DO 274 ! 275 ENDIF 276 ! 277 !== Ice age ==! 278 IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN 186 279 ! clem: in theory we should use the formulation below to advect the ice age, but the code is unable to deal with 187 280 ! fields that do not depend on volume (here oa_i depends on concentration). It creates abnormal ages that … … 189 282 !!zhvar(:,:,:) = poa_i(:,:,:) * z1_ai(:,:,:) 190 283 !!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 age284 ! set u*A=u for advection of ice age 192 285 DO jl = 1, jpl 193 286 zua_ho(:,:,jl) = zudy(:,:) … … 195 288 END DO 196 289 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 ) 290 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho, zva_ho, zcu_box, zcv_box, & 291 & poa_i, poa_i ) 198 292 zamsk = 0._wp 199 293 ENDIF 200 294 ! 201 IF ( ln_pnd_H12 ) THEN !-- melt ponds 202 ! set u*a=u for advection of Ap only 295 !== melt ponds ==! 296 IF ( ln_pnd_H12 ) THEN 297 ! set u*A=u for advection of Ap only 203 298 DO jl = 1, jpl 204 299 zua_ho(:,:,jl) = zudy(:,:) 205 300 zva_ho(:,:,jl) = zvdx(:,:) 206 301 END DO 207 ! 302 ! fraction 208 303 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 304 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 305 & pa_ip, pa_ip, zua_ups, zva_ups, zua_ho , zva_ho ) 210 306 zamsk = 0._wp 211 ! 307 ! volume 212 308 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 309 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 310 & zhvar, pv_ip, zua_ups, zva_ups ) 214 311 ENDIF 215 312 ! 313 !== Open water area ==! 216 314 zati2(:,:) = SUM( pa_i(:,:,:), dim=3 ) 217 315 DO jj = 2, jpjm1 218 316 DO ji = fs_2, fs_jpim1 219 pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & !-- Open water area317 pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & 220 318 & - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 221 319 END DO … … 223 321 CALL lbc_lnk( 'icedyn_adv_umx', pato_i(:,:), 'T', 1. ) 224 322 ! 323 ! 324 ! --- Ensure non-negative fields and in-bound thicknesses --- ! 325 ! 326 ! Remove negative values (conservation is ensured) 327 ! (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 328 CALL ice_var_zapneg( pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 329 ! 330 ! Make sure ice thickness is not too big 331 ! (because ice thickness can be too large where ice concentration is very small) 332 CALL Hbig( zhi_max, zhs_max, zhip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 333 225 334 END DO 226 335 ! … … 228 337 229 338 230 SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pua_ho, pva_ho ) 339 SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, & 340 & pt, ptc, pua_ups, pva_ups, pua_ho, pva_ho ) 231 341 !!---------------------------------------------------------------------- 232 342 !! *** ROUTINE adv_umx *** … … 235 345 !! tracers and add it to the general trend of tracer equations 236 346 !! 237 !! ** Method : - calculate upstream fluxes and upstream solution for tracer H347 !! ** Method : - calculate upstream fluxes and upstream solution for tracers V/A(=H) etc 238 348 !! - calculate tracer H at u and v points (Ultimate) 239 !! - calculate the high order fluxes using alterning directions (Macho ?)349 !! - calculate the high order fluxes using alterning directions (Macho) 240 350 !! - 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 351 !! - convert this tracer flux to a "volume" flux (uH -> uV) 352 !! - apply a limiter a second time on the volumes fluxes (nonosc_ice) 353 !! - calculate the high order solution for V 243 354 !! 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) 355 !! ** Action : solve 3 equations => a) dA/dt = -div(uA) 356 !! b) dV/dt = -div(uV) using dH/dt = -u.grad(H) 357 !! c) dVS/dt = -div(uVS) using either dHS/dt = -u.grad(HS) or dS/dt = -u.grad(S) 250 358 !! 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). 359 !! in eq. b), - fluxes uH are evaluated (with UMx) and limited with nonosc_ice. This step is necessary to get a good H. 360 !! - then we convert this flux to a "volume" flux this way => uH * uA / u 361 !! where uA is the flux from eq. a) 362 !! this "volume" flux is also limited with nonosc_ice (otherwise overshoots can occur) 363 !! - at last we estimate dV/dt = -div(uH * uA / u) 364 !! 365 !! in eq. c), one can solve the equation for S (ln_advS=T), then dVS/dt = -div(uV * uS / u) 366 !! or for HS (ln_advS=F), then dVS/dt = -div(uA * uHS / u) 367 !! 368 !! ** Note : - this method can lead to tiny negative V (-1.e-20) => set it to 0 while conserving mass etc. 369 !! - At the ice edge, Ultimate scheme can lead to: 370 !! 1) negative interpolated tracers at u-v points 371 !! 2) non-zero interpolated tracers at u-v points eventhough there is no ice and velocity is outward 372 !! Solution for 1): apply an upstream scheme when it occurs. A better solution would be to degrade the order of 373 !! the scheme automatically by applying a mask of the ice cover inside Ultimate (not done). 374 !! Solution for 2): we set it to 0 in this case 255 375 !! - Eventhough 1D tests give very good results (typically the one from Schar & Smolarkiewiecz), the 2D is less good. 256 376 !! 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.377 !! work on H (and not V). It is partly related to the multi-category approach 258 378 !! 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 379 !! concentration is small). Since we do not limit S and T, large values can occur at the edge but it does not really matter 380 !! since sv_i and e_i are still good. 381 !!---------------------------------------------------------------------- 382 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) 383 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 384 INTEGER , INTENT(in ) :: jt ! number of sub-iteration 385 INTEGER , INTENT(in ) :: kt ! number of iteration 386 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 387 REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu , pv ! 2 ice velocity components => u*e2 388 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: puc , pvc ! 2 ice velocity components => u*e2 or u*a*e2u 389 REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pubox, pvbox ! upstream velocity 390 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pt ! tracer field 391 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: ptc ! tracer content field 392 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(inout), OPTIONAL :: pua_ups, pva_ups ! upstream u*a fluxes 393 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out), OPTIONAL :: pua_ho, pva_ho ! high order u*a fluxes 273 394 ! 274 395 INTEGER :: ji, jj, jl ! dummy loop indices … … 303 424 DO jj = 1, jpjm1 304 425 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)426 IF( ABS( pu(ji,jj) ) > epsi10 ) THEN 427 zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc (ji,jj,jl) / pu(ji,jj) 428 zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * pua_ups(ji,jj,jl) / pu(ji,jj) 308 429 ELSE 309 430 zfu_ho (ji,jj,jl) = 0._wp … … 311 432 ENDIF 312 433 ! 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)434 IF( ABS( pv(ji,jj) ) > epsi10 ) THEN 435 zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc (ji,jj,jl) / pv(ji,jj) 436 zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pva_ups(ji,jj,jl) / pv(ji,jj) 316 437 ELSE 317 438 zfv_ho (ji,jj,jl) = 0._wp … … 321 442 END DO 322 443 END DO 444 445 ! the new "volume" fluxes must also be "flux corrected" 446 ! thus we calculate the upstream solution and apply a limiter again 447 DO jl = 1, jpl 448 DO jj = 2, jpjm1 449 DO ji = fs_2, fs_jpim1 450 ztra = - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) + zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) 451 ! 452 zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1) 453 END DO 454 END DO 455 END DO 456 CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T', 1. ) 457 ! 458 IF ( kn_limiter == 1 ) THEN 459 CALL nonosc_ice( 1._wp, pdt, pu, pv, ptc, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho ) 460 ELSEIF( kn_limiter == 2 .OR. kn_limiter == 3 ) THEN 461 CALL limiter_x( pdt, pu, ptc, zfu_ups, zfu_ho ) 462 CALL limiter_y( pdt, pv, ptc, zfv_ups, zfv_ho ) 463 ENDIF 464 ! 323 465 ENDIF 324 ! --ho 325 ! in case of advection of A: output u*a 326 ! ------------------------------------- 466 ! --ho --ups 467 ! in case of advection of A: output u*a and u*a 468 ! ----------------------------------------------- 327 469 IF( PRESENT( pua_ho ) ) THEN 328 470 DO jl = 1, jpl 329 471 DO jj = 1, jpjm1 330 472 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 473 pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 474 pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) 475 END DO 334 476 END DO 335 477 END DO … … 609 751 ! 610 752 ! !-- ultimate interpolation of pt at u-point --! 611 CALL ultimate_x( kn_umx, pdt, pt, pu, zt_u, pfu_ho )753 CALL ultimate_x( pamsk, kn_umx, pdt, pt, pu, zt_u, pfu_ho ) 612 754 ! !-- limiter in x --! 613 755 IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) … … 619 761 & + pt (ji,jj,jl) * ( pu (ji,jj ) - pu (ji-1,jj ) ) * r1_e1e2t(ji,jj) & 620 762 & * pamsk & 621 & ) * pdt ) * tmask(ji,jj,1) 763 & ) * pdt ) * tmask(ji,jj,1) 764 !!clem test 765 !!zpt(ji,jj,jl) = MAX( 0._wp, zpt(ji,jj,jl) ) 766 !!clem test 622 767 END DO 623 768 END DO … … 627 772 ! !-- ultimate interpolation of pt at v-point --! 628 773 IF( ll_hoxy ) THEN 629 CALL ultimate_y( kn_umx, pdt, zpt, pv, zt_v, pfv_ho )774 CALL ultimate_y( pamsk, kn_umx, pdt, zpt, pv, zt_v, pfv_ho ) 630 775 ELSE 631 CALL ultimate_y( kn_umx, pdt, pt , pv, zt_v, pfv_ho )776 CALL ultimate_y( pamsk, kn_umx, pdt, pt , pv, zt_v, pfv_ho ) 632 777 ENDIF 633 778 ! !-- limiter in y --! … … 638 783 ! 639 784 ! !-- ultimate interpolation of pt at v-point --! 640 CALL ultimate_y( kn_umx, pdt, pt, pv, zt_v, pfv_ho )785 CALL ultimate_y( pamsk, kn_umx, pdt, pt, pv, zt_v, pfv_ho ) 641 786 ! !-- limiter in y --! 642 787 IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) … … 649 794 & * pamsk & 650 795 & ) * pdt ) * tmask(ji,jj,1) 796 !!clem test 797 !!zpt(ji,jj,jl) = MAX( 0._wp, zpt(ji,jj,jl) ) 798 !!clem test 651 799 END DO 652 800 END DO … … 656 804 ! !-- ultimate interpolation of pt at u-point --! 657 805 IF( ll_hoxy ) THEN 658 CALL ultimate_x( kn_umx, pdt, zpt, pu, zt_u, pfu_ho )806 CALL ultimate_x( pamsk, kn_umx, pdt, zpt, pu, zt_u, pfu_ho ) 659 807 ELSE 660 CALL ultimate_x( kn_umx, pdt, pt , pu, zt_u, pfu_ho )808 CALL ultimate_x( pamsk, kn_umx, pdt, pt , pu, zt_u, pfu_ho ) 661 809 ENDIF 662 810 ! !-- limiter in x --! … … 670 818 671 819 672 SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, pt_u, pfu_ho )820 SUBROUTINE ultimate_x( pamsk, kn_umx, pdt, pt, pu, pt_u, pfu_ho ) 673 821 !!--------------------------------------------------------------------- 674 822 !! *** ROUTINE ultimate_x *** … … 680 828 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 681 829 !!---------------------------------------------------------------------- 830 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) 682 831 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 683 832 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step … … 688 837 ! 689 838 INTEGER :: ji, jj, jl ! dummy loop indices 690 REAL(wp) :: zcu, zdx2, zdx4 ! - -839 REAL(wp) :: zcu, zdx2, zdx4, zvi_cen2 ! - - 691 840 REAL(wp), DIMENSION(jpi,jpj,jpl) :: ztu1, ztu2, ztu3, ztu4 692 841 !!---------------------------------------------------------------------- … … 799 948 END SELECT 800 949 ! 950 ! if there is an outward velocity in a grid cell where there is no ice initially (typically at the ice edge), 951 ! interpolated T at u/v points can be non-zero while it should 952 ! (because of the high order of the advection scheme). Thus set it to 0 in this case 953 IF( ll_icedge ) THEN 954 DO jl = 1, jpl 955 DO jj = 1, jpjm1 956 DO ji = 1, fs_jpim1 957 IF( pt(ji,jj,jl) <= 0._wp .AND. pu(ji,jj) >= 0._wp ) THEN 958 pt_u(ji,jj,jl) = 0._wp 959 ENDIF 960 END DO 961 END DO 962 END DO 963 ENDIF 964 ! 801 965 ! if pt at u-point is negative then use the upstream value 802 966 ! this should not be necessary if a proper sea-ice mask is set in Ultimate … … 806 970 DO jj = 1, jpjm1 807 971 DO ji = 1, fs_jpim1 808 IF( pt_u(ji,jj,jl) < 0._wp ) THEN 972 zvi_cen2 = 0.5_wp * ( v_i(ji+1,jj,jl) + v_i(ji,jj,jl) ) 973 IF( pt_u(ji,jj,jl) < 0._wp .OR. ( zvi_cen2 < epsi06 .AND. pamsk == 0._wp ) ) THEN 809 974 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & 810 975 & - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) … … 826 991 827 992 828 SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pt_v, pfv_ho )993 SUBROUTINE ultimate_y( pamsk, kn_umx, pdt, pt, pv, pt_v, pfv_ho ) 829 994 !!--------------------------------------------------------------------- 830 995 !! *** ROUTINE ultimate_y *** … … 836 1001 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 837 1002 !!---------------------------------------------------------------------- 1003 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) 838 1004 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 839 1005 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step … … 844 1010 ! 845 1011 INTEGER :: ji, jj, jl ! dummy loop indices 846 REAL(wp) :: zcv, zdy2, zdy4 ! - -1012 REAL(wp) :: zcv, zdy2, zdy4, zvi_cen2 ! - - 847 1013 REAL(wp), DIMENSION(jpi,jpj,jpl) :: ztv1, ztv2, ztv3, ztv4 848 1014 !!---------------------------------------------------------------------- … … 952 1118 END SELECT 953 1119 ! 1120 ! if there is an outward velocity in a grid cell where there is no ice initially (typically at the ice edge), 1121 ! interpolated T at u/v points can be non-zero while it should 1122 ! (because of the high order of the advection scheme). Thus set it to 0 in this case 1123 IF( ll_icedge ) THEN 1124 DO jl = 1, jpl 1125 DO jj = 1, jpjm1 1126 DO ji = 1, fs_jpim1 1127 IF( pt(ji,jj,jl) <= 0._wp .AND. pv(ji,jj) >= 0._wp ) THEN 1128 pt_v(ji,jj,jl) = 0._wp 1129 ENDIF 1130 END DO 1131 END DO 1132 END DO 1133 ENDIF 1134 ! 954 1135 ! if pt at v-point is negative then use the upstream value 955 1136 ! this should not be necessary if a proper sea-ice mask is set in Ultimate … … 959 1140 DO jj = 1, jpjm1 960 1141 DO ji = 1, fs_jpim1 961 IF( pt_v(ji,jj,jl) < 0._wp ) THEN 1142 zvi_cen2 = 0.5_wp * ( v_i(ji,jj+1,jl) + v_i(ji,jj,jl) ) 1143 IF( pt_v(ji,jj,jl) < 0._wp .OR. ( zvi_cen2 < epsi06 .AND. pamsk == 0._wp ) ) THEN 962 1144 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) & 963 1145 & - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) … … 1102 1284 ! 1103 1285 ! ! 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 1286 ! 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 !!!) 1287 IF( zpos > epsi10 ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt 1288 ELSE ; zbetup(ji,jj,jl) = 0._wp ! zbig 1106 1289 ENDIF 1107 1290 ! 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 ! zbig1291 IF( zneg > epsi10 ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt 1292 ELSE ; zbetdo(ji,jj,jl) = 0._wp ! zbig 1110 1293 ENDIF 1111 1294 ! … … 1148 1331 END DO 1149 1332 END DO 1150 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 1333 1165 1334 END DO … … 1358 1527 END SUBROUTINE limiter_y 1359 1528 1529 1530 SUBROUTINE Hbig( phi_max, phs_max, phip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 1531 !!------------------------------------------------------------------- 1532 !! *** ROUTINE Hbig *** 1533 !! 1534 !! ** Purpose : Thickness correction in case advection scheme creates 1535 !! abnormally tick ice or snow 1536 !! 1537 !! ** Method : 1- check whether ice thickness is larger than the surrounding 9-points 1538 !! (before advection) and reduce it by adapting ice concentration 1539 !! 2- check whether snow thickness is larger than the surrounding 9-points 1540 !! (before advection) and reduce it by sending the excess in the ocean 1541 !! 3- check whether snow load deplets the snow-ice interface below sea level$ 1542 !! and reduce it by sending the excess in the ocean 1543 !! 4- correct pond fraction to avoid a_ip > a_i 1544 !! 1545 !! ** input : Max thickness of the surrounding 9-points 1546 !!------------------------------------------------------------------- 1547 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: phi_max, phs_max, phip_max ! max ice thick from surrounding 9-pts 1548 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip 1549 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s 1550 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i 1551 ! 1552 INTEGER :: ji, jj, jk, jl ! dummy loop indices 1553 REAL(wp) :: zhip, zhi, zhs, zvs_excess, zfra 1554 REAL(wp), DIMENSION(jpi,jpj) :: zswitch 1555 !!------------------------------------------------------------------- 1556 ! 1557 ! 1558 DO jl = 1, jpl 1559 1560 DO jj = 1, jpj 1561 DO ji = 1, jpi 1562 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 1563 ! 1564 ! ! -- check h_ip -- ! 1565 ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 1566 IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 1567 zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 1568 IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 1569 pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl) 1570 ENDIF 1571 ENDIF 1572 ! 1573 ! ! -- check h_i -- ! 1574 ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i 1575 zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl) 1576 IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 1577 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) 1578 ENDIF 1579 ! 1580 ! ! -- check h_s -- ! 1581 ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean 1582 zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl) 1583 IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 1584 zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 ) 1585 ! 1586 wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * r1_rdtice 1587 hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0 1588 ! 1589 pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 1590 pv_s(ji,jj,jl) = pa_i(ji,jj,jl) * phs_max(ji,jj,jl) 1591 ENDIF 1592 ! 1593 ! ! -- check snow load -- ! 1594 ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean 1595 ! this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin) 1596 ! this imposed mini can artificially make the snow very thick (if concentration decreases drastically) 1597 zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 1598 IF( zvs_excess > 0._wp ) THEN 1599 zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 ) 1600 wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * r1_rdtice 1601 hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0 1602 ! 1603 pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 1604 pv_s(ji,jj,jl) = pv_s(ji,jj,jl) - zvs_excess 1605 ENDIF 1606 1607 ENDIF 1608 END DO 1609 END DO 1610 END DO 1611 ! !-- correct pond fraction to avoid a_ip > a_i 1612 WHERE( pa_ip(:,:,:) > pa_i(:,:,:) ) pa_ip(:,:,:) = pa_i(:,:,:) 1613 ! 1614 ! 1615 END SUBROUTINE Hbig 1616 1360 1617 #else 1361 1618 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
