| [6] | 1 | SUBROUTINE ice_sal_diff(nlay_i,kideb,kiut) |
|---|
| 2 | |
|---|
| 3 | !!------------------------------------------------------------------ |
|---|
| 4 | !! *** ROUTINE ice_sal_diff *** |
|---|
| 5 | !! |
|---|
| 6 | !! ** Purpose : |
|---|
| 7 | !! This routine computes new salinities in the ice |
|---|
| 8 | !! |
|---|
| 9 | !! ** Method : Vertical salinity profile computation |
|---|
| 10 | !! Resolves brine transport equation |
|---|
| 11 | !! |
|---|
| 12 | !! ** Steps |
|---|
| 13 | !! |
|---|
| 14 | !! ** Arguments |
|---|
| 15 | !! |
|---|
| 16 | !! ** Inputs / Outputs |
|---|
| 17 | !! |
|---|
| 18 | !! ** External |
|---|
| 19 | !! |
|---|
| 20 | !! ** References : Vancop. et al., 2008 |
|---|
| 21 | !! |
|---|
| 22 | !! ** History : |
|---|
| 23 | !! (06-2003) Martin Vancop. LIM1D |
|---|
| 24 | !! (06-2008) Martin Vancop. BIO-LIM |
|---|
| 25 | !! (09-2008) Martin Vancop. Explicit gravity drainage |
|---|
| 26 | !! |
|---|
| 27 | !!------------------------------------------------------------------ |
|---|
| 28 | |
|---|
| 29 | USE lib_fortran |
|---|
| 30 | |
|---|
| 31 | INCLUDE 'type.com' |
|---|
| 32 | INCLUDE 'para.com' |
|---|
| 33 | INCLUDE 'const.com' |
|---|
| 34 | INCLUDE 'ice.com' |
|---|
| 35 | INCLUDE 'thermo.com' |
|---|
| 36 | |
|---|
| 37 | REAL(8), DIMENSION(nlay_i) :: |
|---|
| 38 | & z_ms_i , !: mass of salt times thickness |
|---|
| 39 | & z_sbr_i !: brine salinity |
|---|
| 40 | |
|---|
| 41 | REAL(8), DIMENSION(nlay_i) :: !: dummy factors for tracer equation |
|---|
| 42 | & za , !: all |
|---|
| 43 | & zb , !: gravity drainage |
|---|
| 44 | & zc , !: upward advective flow |
|---|
| 45 | & ze , !: downward advective flow |
|---|
| 46 | & zind , !: independent term in the tridiag system |
|---|
| 47 | & zindtbis , !: |
|---|
| 48 | & zdiagbis !: |
|---|
| 49 | |
|---|
| 50 | REAL(8), DIMENSION(nlay_i,3) :: !: dummy factors for tracer equation |
|---|
| 51 | & ztrid !: tridiagonal matrix |
|---|
| 52 | |
|---|
| 53 | REAL(8) :: |
|---|
| 54 | & zdummy1 , !: dummy factors |
|---|
| 55 | & zdummy2 , !: |
|---|
| 56 | & zdummy3 , !: |
|---|
| 57 | & zswitchs , !: switch for summer drainage |
|---|
| 58 | & zeps = 1.0e-20 !: numerical limit |
|---|
| 59 | |
|---|
| 60 | ! Rayleigh number computation |
|---|
| 61 | REAL(8) :: |
|---|
| 62 | & ze_i_min , !: minimum brine volume |
|---|
| 63 | & zcp , !: temporary scalar for sea ice specific heat |
|---|
| 64 | & zk , !: temporary scalar for sea ice thermal conductivity |
|---|
| 65 | & zalphara !: multiplicator for diffusivity |
|---|
| 66 | |
|---|
| 67 | REAL(8), DIMENSION(nlay_i) :: |
|---|
| 68 | & zsigma , !: brine salinity at layer interfaces |
|---|
| 69 | & zperm , !: permeability |
|---|
| 70 | & zpermin , !: minimum permeability |
|---|
| 71 | & zrhodiff , !: density difference |
|---|
| 72 | & zlevel , !: height of the water column |
|---|
| 73 | & zthdiff !: thermal diffusivity |
|---|
| 74 | |
|---|
| 75 | INTEGER :: |
|---|
| 76 | & layer2 , !: layer loop index |
|---|
| 77 | & indtr !: index of tridiagonal system |
|---|
| 78 | |
|---|
| 79 | CHARACTER(len=4) :: |
|---|
| 80 | & bc = 'conc' !: Boundary condition 'conc' or 'flux' |
|---|
| 81 | |
|---|
| 82 | REAL(8) :: |
|---|
| 83 | & z_ms_i_ini , !: initial mass of salt |
|---|
| 84 | & z_ms_i_fin , !: final mass of salt |
|---|
| 85 | & z_fs_b , !: basal flux of salt |
|---|
| 86 | & z_fs_su , !: surface flux of salt |
|---|
| 87 | & z_dms_i !: mass variation |
|---|
| 88 | |
|---|
| 89 | LOGICAL :: |
|---|
| 90 | & ln_write , |
|---|
| 91 | & ln_con , |
|---|
| 92 | & ln_sal |
|---|
| 93 | |
|---|
| 94 | ln_write = .TRUE. ! write outputs |
|---|
| 95 | ln_con = .TRUE. ! conservation check |
|---|
| 96 | ln_sal = .TRUE. ! compute salinity variations or not |
|---|
| 97 | |
|---|
| 98 | IF ( ln_write ) THEN |
|---|
| 99 | WRITE(numout,*) |
|---|
| 100 | WRITE(numout,*) ' ** ice_sal_diff : ' |
|---|
| 101 | WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~ ' |
|---|
| 102 | WRITE(numout,*) ' ln_sal = ', ln_sal |
|---|
| 103 | WRITE(numout,*) ' ln_grd = ', ln_grd |
|---|
| 104 | WRITE(numout,*) ' ln_flu = ', ln_flu |
|---|
| 105 | WRITE(numout,*) ' ln_flo = ', ln_flo |
|---|
| 106 | ENDIF |
|---|
| 107 | WRITE(numout,*) " nlay_i : ", nlay_i |
|---|
| 108 | |
|---|
| 109 | IF ( ln_sal ) THEN |
|---|
| 110 | ! |
|---|
| 111 | !------------------------------------------------------------------------------| |
|---|
| 112 | ! 1) Initialization |
|---|
| 113 | !------------------------------------------------------------------------------| |
|---|
| 114 | ! |
|---|
| 115 | IF ( ln_write ) THEN |
|---|
| 116 | WRITE(numout,*) ' - Initialization ... ' |
|---|
| 117 | ENDIF |
|---|
| 118 | |
|---|
| 119 | DO 10 ji = kideb, kiut |
|---|
| 120 | |
|---|
| 121 | ! brine diffusivity |
|---|
| 122 | diff_br(:) = 0.0 |
|---|
| 123 | |
|---|
| 124 | !--------------------------- |
|---|
| 125 | ! Brine volume and salinity |
|---|
| 126 | !--------------------------- |
|---|
| 127 | DO layer = 1, nlay_i |
|---|
| 128 | e_i_b(layer) = - tmut * s_i_b(ji,layer) / ( t_i_b(ji,layer) |
|---|
| 129 | & - tpw ) |
|---|
| 130 | z_sbr_i(layer) = s_i_b(ji,layer) / e_i_b(layer) |
|---|
| 131 | END DO |
|---|
| 132 | |
|---|
| 133 | !-------------------- |
|---|
| 134 | ! Conservation check |
|---|
| 135 | !-------------------- |
|---|
| 136 | IF ( ln_con ) THEN |
|---|
| 137 | CALL ice_sal_column( kideb , kiut , z_ms_i_ini , |
|---|
| 138 | & s_i_b(1,1:nlay_i), |
|---|
| 139 | & deltaz_i_phy, nlay_i, .FALSE. ) |
|---|
| 140 | ENDIF ! ln_con |
|---|
| 141 | |
|---|
| 142 | IF ( ln_write ) THEN |
|---|
| 143 | WRITE(numout,*) ' nlay_i : ', nlay_i |
|---|
| 144 | WRITE(numout,*) ' kideb : ', kideb |
|---|
| 145 | WRITE(numout,*) ' kiut : ', kiut |
|---|
| 146 | WRITE(numout,*) ' ' |
|---|
| 147 | WRITE(numout,*) ' deltaz_i_phy : ', ( deltaz_i_phy(layer), |
|---|
| 148 | & layer = 1, nlay_i ) |
|---|
| 149 | WRITE(numout,*) ' z_i_phy : ', ( z_i_phy(layer), |
|---|
| 150 | & layer = 1, nlay_i ) |
|---|
| 151 | WRITE(numout,*) ' s_i_b : ', ( s_i_b (ji,layer), |
|---|
| 152 | & layer = 1, nlay_i ) |
|---|
| 153 | WRITE(numout,*) ' t_i_b : ', ( t_i_b (ji,layer), |
|---|
| 154 | & layer = 1, nlay_i ) |
|---|
| 155 | WRITE(numout,*) ' e_i_b : ', ( e_i_b (layer), |
|---|
| 156 | & layer = 1, nlay_i ) |
|---|
| 157 | WRITE(numout,*) ' z_sbr_i : ', ( z_sbr_i (layer), |
|---|
| 158 | & layer = 1, nlay_i ) |
|---|
| 159 | WRITE(numout,*) |
|---|
| 160 | ENDIF ! ln_write |
|---|
| 161 | |
|---|
| 162 | 10 CONTINUE |
|---|
| 163 | |
|---|
| 164 | ! |
|---|
| 165 | !------------------------------------------------------------------------------| |
|---|
| 166 | ! 2) Rayleigh-number-based diffusivity |
|---|
| 167 | !------------------------------------------------------------------------------| |
|---|
| 168 | ! |
|---|
| 169 | ! Diffusivity is a function of the local Rayleigh number |
|---|
| 170 | ! see Notz and Worster, JGR 2008 |
|---|
| 171 | ! Diffusivity, layer represents the interface' |
|---|
| 172 | ! between layer and layer-1 ' |
|---|
| 173 | ! |
|---|
| 174 | IF ( ln_write ) THEN |
|---|
| 175 | WRITE(numout,*) ' - Rayleigh-number based diffusivity ... ' |
|---|
| 176 | WRITE(numout,*) ' ' |
|---|
| 177 | ENDIF |
|---|
| 178 | |
|---|
| 179 | DO 20 ji = kideb, kiut |
|---|
| 180 | |
|---|
| 181 | !----------------------------------------- |
|---|
| 182 | ! Brine salinity between layer interfaces |
|---|
| 183 | !----------------------------------------- |
|---|
| 184 | DO layer = 1, nlay_i - 1 |
|---|
| 185 | zdummy1 = t_i_b(ji,layer) + deltaz_i_phy(layer) / 2. * |
|---|
| 186 | & ( t_i_b(ji,layer+1) - t_i_b(ji,layer) ) / |
|---|
| 187 | & ( z_i_phy(layer+1) - z_i_phy(layer) ) - tpw |
|---|
| 188 | zsigma(layer) = - zdummy1 / tmut |
|---|
| 189 | END DO |
|---|
| 190 | zsigma(nlay_i) = - ( t_i_b(ji,nlay_i) - tpw ) / tmut |
|---|
| 191 | |
|---|
| 192 | !-------------------- |
|---|
| 193 | ! Density difference |
|---|
| 194 | !-------------------- |
|---|
| 195 | DO layer = 1, nlay_i |
|---|
| 196 | zrhodiff(layer) = - beta_ocs * ( oce_sal - zsigma(layer) ) |
|---|
| 197 | END DO |
|---|
| 198 | |
|---|
| 199 | !------------------------------------------ |
|---|
| 200 | ! Minimum permeability under current level |
|---|
| 201 | !------------------------------------------ |
|---|
| 202 | DO layer = 1, nlay_i |
|---|
| 203 | ze_i_min = 99999.0 |
|---|
| 204 | DO layer2 = layer, nlay_i |
|---|
| 205 | ze_i_min = MIN( ze_i_min , e_i_b(layer2) ) |
|---|
| 206 | zpermin(layer) = 1.0e-17 * ( ( 1000. * ze_i_min )**3.1 ) |
|---|
| 207 | END DO |
|---|
| 208 | END DO ! layer |
|---|
| 209 | |
|---|
| 210 | !------------------------------------------------ |
|---|
| 211 | ! length of the water column under current level |
|---|
| 212 | !------------------------------------------------ |
|---|
| 213 | DO layer = nlay_i, 1, -1 |
|---|
| 214 | zlevel(layer) = 0.0 |
|---|
| 215 | DO layer2 = layer, nlay_i |
|---|
| 216 | zlevel(layer) = zlevel(layer) + deltaz_i_phy(layer2) |
|---|
| 217 | END DO |
|---|
| 218 | END DO |
|---|
| 219 | zlevel(nlay_i) = deltaz_i_phy(nlay_i) / 2.0 |
|---|
| 220 | |
|---|
| 221 | !--------------------- |
|---|
| 222 | ! Thermal diffusivity |
|---|
| 223 | !--------------------- |
|---|
| 224 | zkimin = 0.1 |
|---|
| 225 | DO layer = 1, nlay_i - 1 |
|---|
| 226 | zdummy1 = t_i_b(ji,layer) + deltaz_i_phy(layer) / 2. * |
|---|
| 227 | & ( t_i_b(ji,layer+1) - t_i_b(ji,layer) ) / |
|---|
| 228 | & ( z_i_phy(layer+1) - z_i_phy(layer) ) - tpw |
|---|
| 229 | zdummy2 = s_i_b(ji,layer) + deltaz_i_phy(layer) / 2. * |
|---|
| 230 | & ( s_i_b(ji,layer+1) - s_i_b(ji,layer) ) / |
|---|
| 231 | & ( z_i_phy(layer+1) - z_i_phy(layer) ) |
|---|
| 232 | zcp = cpg + lfus * tmut * zdummy2 / |
|---|
| 233 | & MAX( zdummy1 * zdummy1 , zeps ) |
|---|
| 234 | zk = xkg + betak1 * zdummy2 / |
|---|
| 235 | & MIN( -zeps , zdummy1 ) - betak2 * zdummy1 |
|---|
| 236 | zk = MAX( zk, zkimin ) |
|---|
| 237 | zthdiff(layer) = zk / ( rhog * zcp ) |
|---|
| 238 | END DO |
|---|
| 239 | |
|---|
| 240 | zcp = cpg + lfus * tmut * s_i_b(ji,nlay_i) / |
|---|
| 241 | & MAX( ( t_i_b(ji,nlay_i) - tpw ) * ( t_i_b(ji,nlay_i) - tpw ), |
|---|
| 242 | & zeps ) |
|---|
| 243 | zk = xkg + betak1 * s_i_b(ji,nlay_i) / |
|---|
| 244 | & MIN( -zeps , t_i_b(ji,nlay_i) - tpw ) |
|---|
| 245 | & - betak2 * ( t_i_b(ji,nlay_i) - tpw ) |
|---|
| 246 | zk = MAX( zk, zkimin ) |
|---|
| 247 | zthdiff(nlay_i) = zk / ( rhog * zcp ) |
|---|
| 248 | |
|---|
| 249 | !----------------- |
|---|
| 250 | ! Rayleigh number |
|---|
| 251 | !----------------- |
|---|
| 252 | DO layer = 1, nlay_i |
|---|
| 253 | rayleigh(layer) = gpes * MAX(zrhodiff(layer),0.0) * |
|---|
| 254 | & zpermin(layer) * zlevel(layer) / |
|---|
| 255 | & ( zthdiff(layer) * visc_br ) |
|---|
| 256 | END DO |
|---|
| 257 | |
|---|
| 258 | !------------------- |
|---|
| 259 | ! Brine Diffusivity |
|---|
| 260 | !------------------- |
|---|
| 261 | DO layer = 1, nlay_i |
|---|
| 262 | zalphara = ( TANH( ra_smooth * ( rayleigh(layer) - ra_c ) ) |
|---|
| 263 | & + 1 ) / 2.0 |
|---|
| 264 | diff_br(layer) = ( 1.0 - zalphara ) * d_br_mol + |
|---|
| 265 | & zalphara * ( d_br_tur ) |
|---|
| 266 | IF ( .NOT. ln_grd ) diff_br(layer) = 0. |
|---|
| 267 | END DO |
|---|
| 268 | |
|---|
| 269 | |
|---|
| 270 | IF ( ln_write ) THEN |
|---|
| 271 | WRITE(numout,*) ' zsigma : ', ( zsigma(layer), |
|---|
| 272 | & layer = 1, nlay_i) |
|---|
| 273 | WRITE(numout,*) ' zrhodiff : ', ( zrhodiff(layer), |
|---|
| 274 | & layer = 1, nlay_i ) |
|---|
| 275 | WRITE(numout,*) ' zpermin : ', ( zpermin(layer), |
|---|
| 276 | & layer = 1, nlay_i ) |
|---|
| 277 | WRITE(numout,*) ' zthdiff : ', ( zthdiff(layer), |
|---|
| 278 | & layer = 1, nlay_i ) |
|---|
| 279 | WRITE(numout,*) ' zlevel : ', ( zlevel(layer), |
|---|
| 280 | & layer = 1, nlay_i ) |
|---|
| 281 | WRITE(numout,*) ' rayleigh : ', ( rayleigh(layer), |
|---|
| 282 | & layer = 1, nlay_i ) |
|---|
| 283 | WRITE(numout,*) ' diff_br : ', ( diff_br(layer), |
|---|
| 284 | & layer = 1, nlay_i ) |
|---|
| 285 | WRITE(numout,*) |
|---|
| 286 | ENDIF |
|---|
| 287 | |
|---|
| 288 | 20 CONTINUE |
|---|
| 289 | |
|---|
| 290 | ! |
|---|
| 291 | !------------------------------------------------------------------------------| |
|---|
| 292 | ! 3) Flooding and flushing velocities |
|---|
| 293 | !------------------------------------------------------------------------------| |
|---|
| 294 | ! |
|---|
| 295 | IF ( ln_write ) THEN |
|---|
| 296 | WRITE(numout,*) ' - Flooding and flushing velocities ' |
|---|
| 297 | WRITE(numout,*) ' ' |
|---|
| 298 | ENDIF |
|---|
| 299 | |
|---|
| 300 | DO 30 ji = kideb, kiut |
|---|
| 301 | !----------------------- |
|---|
| 302 | ! Permeability switches |
|---|
| 303 | !----------------------- |
|---|
| 304 | ! Permeability switch = 1 if brine volume fraction > e_thr_flu |
|---|
| 305 | zswitch_per = 1.0 |
|---|
| 306 | zbvmin = 1.0 |
|---|
| 307 | DO layer = 1, nlay_i |
|---|
| 308 | zbvmin = MIN( e_i_b(layer) , zbvmin ) ! minimum brine volume |
|---|
| 309 | END DO |
|---|
| 310 | IF ( zbvmin .LT. e_thr_flu ) zswitch_per = 0.0 |
|---|
| 311 | |
|---|
| 312 | ! Summer switch = 1 if Tsu ge tpw and min brine volume superior than e_thr_flu |
|---|
| 313 | zswitchs = MAX( 0.0, SIGN ( 1.0d0, t_su_b(ji) - tpw ) ) ! 0 if winter, 1 if summer |
|---|
| 314 | zswitchs = zswitchs * zswitch_per |
|---|
| 315 | |
|---|
| 316 | !------------------- |
|---|
| 317 | ! Flooding velocity |
|---|
| 318 | !------------------- |
|---|
| 319 | zdhitot = dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) |
|---|
| 320 | zdhstot = dh_s_melt(ji) |
|---|
| 321 | w_flood = ( ( rho0 - rhog ) / rho0 * zdhitot - |
|---|
| 322 | & rhon / rho0 * zdhstot ) / ddtb * zswitch_per |
|---|
| 323 | IF ( w_flood .GT. 0 ) THEN |
|---|
| 324 | z_flood = 0.0 |
|---|
| 325 | ELSE |
|---|
| 326 | z_flood = 1.0 |
|---|
| 327 | ENDIF |
|---|
| 328 | IF ( .NOT. ln_flo ) w_flood = 0. |
|---|
| 329 | |
|---|
| 330 | !------------------ |
|---|
| 331 | ! Percolating flux |
|---|
| 332 | !------------------ |
|---|
| 333 | ! Percolating flow ( rho dh * beta * switch / rho0 ) |
|---|
| 334 | qsummer = ( - rhog * MIN ( dh_i_surf(ji) , 0.0 ) |
|---|
| 335 | & - rhon * MIN ( dh_s_tot(ji) , 0.0 ) ) |
|---|
| 336 | qsummer = qsummer * flu_beta * zswitchs / 1000.0 ! 1000 is a ref density for brine |
|---|
| 337 | |
|---|
| 338 | w_flush = qsummer / ddtb / e_i_b(1) |
|---|
| 339 | |
|---|
| 340 | IF ( .NOT. ln_flu ) THEN |
|---|
| 341 | w_flush = 0. |
|---|
| 342 | qsummer = 0. |
|---|
| 343 | ENDIF |
|---|
| 344 | |
|---|
| 345 | IF ( ln_write ) THEN |
|---|
| 346 | WRITE(numout,*) ' zswitchs : ', zswitchs |
|---|
| 347 | WRITE(numout,*) ' w_flush : ', w_flush |
|---|
| 348 | WRITE(numout,*) ' w_flood : ', w_flood |
|---|
| 349 | ENDIF |
|---|
| 350 | |
|---|
| 351 | 30 CONTINUE |
|---|
| 352 | ! |
|---|
| 353 | !------------------------------------------------------------------------------| |
|---|
| 354 | ! 4) Compute dummy factors for tracer diffusion equation |
|---|
| 355 | !------------------------------------------------------------------------------| |
|---|
| 356 | ! |
|---|
| 357 | IF ( ln_write ) THEN |
|---|
| 358 | WRITE(numout,*) ' - Compute dummy factors for tracer diffusion' |
|---|
| 359 | WRITE(numout,*) ' ' |
|---|
| 360 | ENDIF |
|---|
| 361 | |
|---|
| 362 | DO 40 ji = kideb, kiut |
|---|
| 363 | !---------------- |
|---|
| 364 | ! za factors |
|---|
| 365 | !---------------- |
|---|
| 366 | DO layer = 1, nlay_i |
|---|
| 367 | za(layer) = ddtb / ( deltaz_i_phy(layer) * e_i_b(layer) ) |
|---|
| 368 | END DO |
|---|
| 369 | |
|---|
| 370 | !-------------------- |
|---|
| 371 | ! zb, zc, ze factors |
|---|
| 372 | !-------------------- |
|---|
| 373 | DO layer = 1, nlay_i - 1 |
|---|
| 374 | ! interpolate brine volume at the interface between layers |
|---|
| 375 | zdummy1 = ( e_i_b(layer + 1 ) - e_i_b(layer) ) / |
|---|
| 376 | & ( z_i_phy(layer + 1) - z_i_phy(layer) ) |
|---|
| 377 | zdummy2 = deltaz_i_phy(layer) / 2.0 |
|---|
| 378 | zdummy3 = e_i_b(layer) + zdummy1 * zdummy2 |
|---|
| 379 | zb(layer) = zdummy3 * diff_br(layer) / |
|---|
| 380 | & ( z_i_phy(layer + 1) - z_i_phy(layer) ) |
|---|
| 381 | zc(layer) = w_flood * zdummy3 * z_flood |
|---|
| 382 | ze(layer) = ( w_flood * ( 1. - z_flood ) + w_flush ) * zdummy3 |
|---|
| 383 | ! old qsummer scheme |
|---|
| 384 | ! ze(layer) = ( w_flood * ( 1. - z_flood ) ) * zdummy3 + |
|---|
| 385 | ! & zswitchs * qsummer / ddtb |
|---|
| 386 | END DO |
|---|
| 387 | |
|---|
| 388 | ! Fixed boundary condition (imposed cc.) |
|---|
| 389 | zb(nlay_i) = 2. * e_i_b(nlay_i) / |
|---|
| 390 | & deltaz_i_phy(nlay_i) * diff_br(nlay_i) |
|---|
| 391 | zc(nlay_i) = w_flood * e_i_b(nlay_i) * z_flood |
|---|
| 392 | ze(nlay_i) = ( w_flood * ( 1. - z_flood ) + w_flush ) * |
|---|
| 393 | & e_i_b(nlay_i) |
|---|
| 394 | ! ze(nlay_i) = ( w_flood * ( 1. - z_flood ) ) * e_i_b(nlay_i) + |
|---|
| 395 | ! & zswitchs * qsummer / ddtb |
|---|
| 396 | |
|---|
| 397 | IF ( ln_write ) THEN |
|---|
| 398 | WRITE(numout,*) |
|---|
| 399 | WRITE(numout,*) ' -Dummy factors ' |
|---|
| 400 | WRITE(numout,*) ' za : ', ( za (layer), |
|---|
| 401 | & layer = 1, nlay_i) |
|---|
| 402 | WRITE(numout,*) ' zb : ', ( zb (layer), |
|---|
| 403 | & layer = 1, nlay_i) |
|---|
| 404 | WRITE(numout,*) ' zc : ', ( zc (layer), |
|---|
| 405 | & layer = 1, nlay_i) |
|---|
| 406 | WRITE(numout,*) ' ze : ', ( ze (layer), |
|---|
| 407 | & layer = 1, nlay_i) |
|---|
| 408 | ENDIF |
|---|
| 409 | |
|---|
| 410 | 40 CONTINUE |
|---|
| 411 | ! |
|---|
| 412 | !----------------------------------------------------------------------- |
|---|
| 413 | ! 5) Tridiagonal system terms for tracer diffusion equation |
|---|
| 414 | !----------------------------------------------------------------------- |
|---|
| 415 | ! |
|---|
| 416 | DO 50 ji = kideb, kiut |
|---|
| 417 | |
|---|
| 418 | !---------------- |
|---|
| 419 | ! first equation |
|---|
| 420 | !---------------- |
|---|
| 421 | ztrid(1,1) = 0.0 |
|---|
| 422 | ztrid(1,2) = 1.0 + za(1) * ( zb(1) + ze(1) ) |
|---|
| 423 | ztrid(1,3) = za(1) * ( -zb(1) + zc(1) ) |
|---|
| 424 | zind(1) = z_sbr_i(1) |
|---|
| 425 | |
|---|
| 426 | !----------------- |
|---|
| 427 | ! inner equations |
|---|
| 428 | !----------------- |
|---|
| 429 | DO layer = 2, nlay_i - 1 |
|---|
| 430 | ztrid(layer,1) = - za(layer) * ( zb(layer-1) + ze(layer-1) ) |
|---|
| 431 | ztrid(layer,2) = 1.0 + za(layer) * ( zb(layer) + |
|---|
| 432 | & ze(layer) + zb(layer-1) - zc(layer-1) ) |
|---|
| 433 | ztrid(layer,3) = za(layer) * ( -zb(layer) + zc(layer) ) |
|---|
| 434 | zind(layer) = z_sbr_i(layer) |
|---|
| 435 | END DO |
|---|
| 436 | |
|---|
| 437 | !---------------- |
|---|
| 438 | ! last equation |
|---|
| 439 | !---------------- |
|---|
| 440 | WRITE(numout,*) " nlay_i : ", nlay_i |
|---|
| 441 | ztrid(nlay_i,1) = - za(nlay_i) * ( zb(nlay_i-1) + ze(nlay_i-1) ) |
|---|
| 442 | ztrid(nlay_i,2) = 1.0 + za(nlay_i) * ( zb(nlay_i) + |
|---|
| 443 | & ze(nlay_i) + zb(nlay_i-1) - zc(nlay_i-1) ) |
|---|
| 444 | ztrid(nlay_i,3) = 0. |
|---|
| 445 | zind(nlay_i) = z_sbr_i(nlay_i) + za(nlay_i) * ( zb(nlay_i) |
|---|
| 446 | & - zc(nlay_i) ) * oce_sal |
|---|
| 447 | |
|---|
| 448 | IF ( ln_write ) THEN |
|---|
| 449 | WRITE(numout,*) |
|---|
| 450 | WRITE(numout,*) ' -Tridiag terms, winter ... ' |
|---|
| 451 | WRITE(numout,*) |
|---|
| 452 | DO layer = 1, nlay_i |
|---|
| 453 | WRITE(numout,*) ' layer : ', layer |
|---|
| 454 | WRITE(numout,*) ' ztrid : ', ztrid(layer,1), |
|---|
| 455 | & ztrid(layer,2), ztrid(layer,3) |
|---|
| 456 | WRITE(numout,*) ' zind : ', zind(layer) |
|---|
| 457 | END DO |
|---|
| 458 | ENDIF |
|---|
| 459 | |
|---|
| 460 | 50 CONTINUE |
|---|
| 461 | |
|---|
| 462 | ! |
|---|
| 463 | !----------------------------------------------------------------------- |
|---|
| 464 | ! 6) Solving the tridiagonal system |
|---|
| 465 | !----------------------------------------------------------------------- |
|---|
| 466 | ! |
|---|
| 467 | DO 60 ji = kideb, kiut |
|---|
| 468 | |
|---|
| 469 | ! The tridiagonal system is solved with Gauss elimination |
|---|
| 470 | ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, |
|---|
| 471 | ! McGraw-Hill 1984. |
|---|
| 472 | |
|---|
| 473 | zindtbis(1) = zind(1) |
|---|
| 474 | zdiagbis(1) = ztrid(1,2) |
|---|
| 475 | |
|---|
| 476 | DO layer = 2, nlay_i |
|---|
| 477 | zdiagbis(layer) = ztrid(layer,2) - ztrid(layer,1) * |
|---|
| 478 | & ztrid(layer-1,3) / zdiagbis(layer-1) |
|---|
| 479 | zindtbis(layer) = zind(layer) - ztrid(layer,1) * |
|---|
| 480 | & zindtbis(layer-1) / zdiagbis(layer-1) |
|---|
| 481 | END DO |
|---|
| 482 | |
|---|
| 483 | ! Recover brine salinity |
|---|
| 484 | z_sbr_i(nlay_i) = zindtbis(nlay_i) / zdiagbis(nlay_i) |
|---|
| 485 | DO layer = nlay_i - 1 , 1 , -1 |
|---|
| 486 | z_sbr_i(layer) = ( zindtbis(layer) - ztrid(layer,3)* |
|---|
| 487 | & z_sbr_i(layer+1)) / zdiagbis(layer) |
|---|
| 488 | END DO |
|---|
| 489 | ! Recover ice salinity |
|---|
| 490 | DO layer = 1, nlay_i |
|---|
| 491 | sn_i_b(layer) = z_sbr_i(layer) * e_i_b(layer) |
|---|
| 492 | ! s_i_b(ji,layer) = z_sbr_i(layer) * e_i_b(layer) |
|---|
| 493 | END DO |
|---|
| 494 | |
|---|
| 495 | IF ( ln_write ) THEN |
|---|
| 496 | WRITE(numout,*) |
|---|
| 497 | WRITE(numout,*) ' -Solving the tridiagonal system ... ' |
|---|
| 498 | WRITE(numout,*) |
|---|
| 499 | WRITE(numout,*) ' zdiagbis: ', ( zdiagbis(layer) , |
|---|
| 500 | & layer = 1, nlay_i ) |
|---|
| 501 | WRITE(numout,*) ' zindtbis: ', ( zdiagbis(layer) , |
|---|
| 502 | & layer = 1, nlay_i ) |
|---|
| 503 | WRITE(numout,*) ' z_sbr_i : ', ( z_sbr_i(layer) , |
|---|
| 504 | & layer = 1, nlay_i ) |
|---|
| 505 | WRITE(numout,*) ' sn_i_b : ', ( sn_i_b(layer) , |
|---|
| 506 | & layer = 1, nlay_i ) |
|---|
| 507 | ENDIF |
|---|
| 508 | |
|---|
| 509 | 60 CONTINUE |
|---|
| 510 | ! |
|---|
| 511 | !----------------------------------------------------------------------- |
|---|
| 512 | ! 7) Mass of salt conserved ? |
|---|
| 513 | !----------------------------------------------------------------------- |
|---|
| 514 | ! |
|---|
| 515 | DO 70 ji = kideb, kiut |
|---|
| 516 | |
|---|
| 517 | ! Final mass of salt |
|---|
| 518 | CALL ice_sal_column( kideb , kiut , z_ms_i_fin , |
|---|
| 519 | & sn_i_b(1:nlay_i), |
|---|
| 520 | & deltaz_i_phy, nlay_i, .FALSE. ) |
|---|
| 521 | |
|---|
| 522 | ! Bottom flux ( positive upwards for conservation routine ) |
|---|
| 523 | zfb = - e_i_b(nlay_i) * |
|---|
| 524 | & ( diff_br(nlay_i) * 2.0 / deltaz_i_phy(nlay_i) * |
|---|
| 525 | & ( z_sbr_i(nlay_i) - oce_sal ) + w_flood * ( z_flood * |
|---|
| 526 | & oce_sal + ( 1. - z_flood ) * z_sbr_i(nlay_i) ) ) |
|---|
| 527 | & - e_i_b(nlay_i) * w_flush * z_sbr_i(nlay_i) |
|---|
| 528 | ! & - qsummer * z_sbr_i(nlay_i) / ddtb |
|---|
| 529 | |
|---|
| 530 | fsb = - zfb * rhog / 1000. ! ice-ocean salt flux is positive downwards |
|---|
| 531 | IF ( ln_write ) THEN |
|---|
| 532 | WRITE(numout,*) ' fsb : ', fsb |
|---|
| 533 | WRITE(numout,*) |
|---|
| 534 | ENDIF |
|---|
| 535 | |
|---|
| 536 | ! Surface flux of salt |
|---|
| 537 | zfsu = 0.0 |
|---|
| 538 | |
|---|
| 539 | ! conservation check |
|---|
| 540 | zerror = 1.0e-15 |
|---|
| 541 | CALL ice_sal_conserv(kideb,kiut,'ice_sal_diff : ',zerror, |
|---|
| 542 | & z_ms_i_ini,z_ms_i_fin, |
|---|
| 543 | & zfb , zfsu , ddtb) |
|---|
| 544 | |
|---|
| 545 | 70 CONTINUE |
|---|
| 546 | |
|---|
| 547 | ENDIF ! ln_sal |
|---|
| 548 | ! |
|---|
| 549 | !------------------------------------------------------------------------------| |
|---|
| 550 | ! End of la sous-routine |
|---|
| 551 | WRITE(numout,*) |
|---|
| 552 | END SUBROUTINE |
|---|
![(please configure the [header_logo] section in trac.ini)](/lim1d/chrome/site/your_project_logo.png){kind=logo}
