Changeset 4870 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
- Timestamp:
- 2014-11-18T17:03:01+01:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r4765 r4870 98 98 INTEGER :: ii, ij ! temporary dummy loop index 99 99 INTEGER :: numeq ! current reference number of equation 100 INTEGER :: layer! vertical dummy loop index100 INTEGER :: jk ! vertical dummy loop index 101 101 INTEGER :: nconv ! number of iterations in iterative procedure 102 102 INTEGER :: minnumeqmin, maxnumeqmax … … 188 188 z_i(:,0) = 0._wp ! vert. coord. of the up. lim. of the 1st ice layer 189 189 190 DO layer= 1, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer191 DO ji = kideb , kiut 192 z_s(ji, layer) = z_s(ji,layer-1) + ht_s_b(ji) / REAL( nlay_s )193 END DO 194 END DO 195 196 DO layer= 1, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer197 DO ji = kideb , kiut 198 z_i(ji, layer) = z_i(ji,layer-1) + ht_i_b(ji) / REAL( nlay_i )190 DO jk = 1, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer 191 DO ji = kideb , kiut 192 z_s(ji,jk) = z_s(ji,jk-1) + ht_s_b(ji) / REAL( nlay_s ) 193 END DO 194 END DO 195 196 DO jk = 1, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer 197 DO ji = kideb , kiut 198 z_i(ji,jk) = z_i(ji,jk-1) + ht_i_b(ji) / REAL( nlay_i ) 199 199 END DO 200 200 END DO … … 249 249 END DO 250 250 251 DO layer= 1, nlay_s ! Radiation through snow251 DO jk = 1, nlay_s ! Radiation through snow 252 252 DO ji = kideb, kiut 253 253 ! ! radiation transmitted below the layer-th snow layer 254 zradtr_s(ji, layer) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,layer) ) ) )254 zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,jk) ) ) ) 255 255 ! ! radiation absorbed by the layer-th snow layer 256 zradab_s(ji, layer) = zradtr_s(ji,layer-1) - zradtr_s(ji,layer)256 zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk) 257 257 END DO 258 258 END DO … … 262 262 END DO 263 263 264 DO layer= 1, nlay_i ! Radiation through ice264 DO jk = 1, nlay_i ! Radiation through ice 265 265 DO ji = kideb, kiut 266 266 ! ! radiation transmitted below the layer-th ice layer 267 zradtr_i(ji, layer) = zradtr_i(ji,0) * EXP( - kappa_i * ( MAX ( 0._wp , z_i(ji,layer) ) ) )267 zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - kappa_i * ( MAX ( 0._wp , z_i(ji,jk) ) ) ) 268 268 ! ! radiation absorbed by the layer-th ice layer 269 zradab_i(ji, layer) = zradtr_i(ji,layer-1) - zradtr_i(ji,layer)269 zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk) 270 270 END DO 271 271 END DO … … 288 288 END DO 289 289 290 DO layer= 1, nlay_s ! Old snow temperature291 DO ji = kideb , kiut 292 ztsold(ji, layer) = t_s_b(ji,layer)293 END DO 294 END DO 295 296 DO layer= 1, nlay_i ! Old ice temperature297 DO ji = kideb , kiut 298 ztiold(ji, layer) = t_i_b(ji,layer)290 DO jk = 1, nlay_s ! Old snow temperature 291 DO ji = kideb , kiut 292 ztsold(ji,jk) = t_s_b(ji,jk) 293 END DO 294 END DO 295 296 DO jk = 1, nlay_i ! Old ice temperature 297 DO ji = kideb , kiut 298 ztiold(ji,jk) = t_i_b(ji,jk) 299 299 END DO 300 300 END DO … … 316 316 ztcond_i(ji,0) = MAX(ztcond_i(ji,0),zkimin) 317 317 END DO 318 DO layer= 1, nlay_i-1318 DO jk = 1, nlay_i-1 319 319 DO ji = kideb , kiut 320 ztcond_i(ji, layer) = rcdic + zbeta*( s_i_b(ji,layer) + s_i_b(ji,layer+1) ) / &321 MIN(-2.0_wp * epsi10, t_i_b(ji, layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)322 ztcond_i(ji, layer) = MAX(ztcond_i(ji,layer),zkimin)320 ztcond_i(ji,jk) = rcdic + zbeta*( s_i_b(ji,jk) + s_i_b(ji,jk+1) ) / & 321 MIN(-2.0_wp * epsi10, t_i_b(ji,jk)+t_i_b(ji,jk+1) - 2.0_wp * rtt) 322 ztcond_i(ji,jk) = MAX(ztcond_i(ji,jk),zkimin) 323 323 END DO 324 324 END DO … … 331 331 ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) 332 332 END DO 333 DO layer= 1, nlay_i-1333 DO jk = 1, nlay_i-1 334 334 DO ji = kideb , kiut 335 ztcond_i(ji, layer) = rcdic + &336 & 0.090_wp * ( s_i_b(ji, layer) + s_i_b(ji,layer+1) ) &337 & / MIN(-2.0_wp * epsi10, t_i_b(ji, layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt) &338 & - 0.0055_wp* ( t_i_b(ji, layer) + t_i_b(ji,layer+1) - 2.0*rtt )339 ztcond_i(ji, layer) = MAX( ztcond_i(ji,layer), zkimin )335 ztcond_i(ji,jk) = rcdic + & 336 & 0.090_wp * ( s_i_b(ji,jk) + s_i_b(ji,jk+1) ) & 337 & / MIN(-2.0_wp * epsi10, t_i_b(ji,jk)+t_i_b(ji,jk+1) - 2.0_wp * rtt) & 338 & - 0.0055_wp* ( t_i_b(ji,jk) + t_i_b(ji,jk+1) - 2.0*rtt ) 339 ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin ) 340 340 END DO 341 341 END DO … … 358 358 END DO 359 359 360 DO layer= 1, nlay_s-1361 DO ji = kideb , kiut 362 zkappa_s(ji, layer) = 2.0 * rcdsn / &360 DO jk = 1, nlay_s-1 361 DO ji = kideb , kiut 362 zkappa_s(ji,jk) = 2.0 * rcdsn / & 363 363 MAX(epsi10,2.0*zh_s(ji)) 364 364 END DO 365 365 END DO 366 366 367 DO layer= 1, nlay_i-1367 DO jk = 1, nlay_i-1 368 368 DO ji = kideb , kiut 369 369 !-- Ice kappa factors 370 zkappa_i(ji, layer) = 2.0*ztcond_i(ji,layer)/ &370 zkappa_i(ji,jk) = 2.0*ztcond_i(ji,jk)/ & 371 371 MAX(epsi10,2.0*zh_i(ji)) 372 372 END DO … … 387 387 !------------------------------------------------------------------------------| 388 388 ! 389 DO layer= 1, nlay_i390 DO ji = kideb , kiut 391 ztitemp(ji, layer) = t_i_b(ji,layer)392 zspeche_i(ji, layer) = cpic + zgamma*s_i_b(ji,layer)/ &393 MAX((t_i_b(ji, layer)-rtt)*(ztiold(ji,layer)-rtt),epsi10)394 zeta_i(ji, layer) = rdt_ice / MAX(rhoic*zspeche_i(ji,layer)*zh_i(ji), &389 DO jk = 1, nlay_i 390 DO ji = kideb , kiut 391 ztitemp(ji,jk) = t_i_b(ji,jk) 392 zspeche_i(ji,jk) = cpic + zgamma*s_i_b(ji,jk)/ & 393 MAX((t_i_b(ji,jk)-rtt)*(ztiold(ji,jk)-rtt),epsi10) 394 zeta_i(ji,jk) = rdt_ice / MAX(rhoic*zspeche_i(ji,jk)*zh_i(ji), & 395 395 epsi10) 396 396 END DO 397 397 END DO 398 398 399 DO layer= 1, nlay_s400 DO ji = kideb , kiut 401 ztstemp(ji, layer) = t_s_b(ji,layer)402 zeta_s(ji, layer) = rdt_ice / MAX(rhosn*cpic*zh_s(ji),epsi10)399 DO jk = 1, nlay_s 400 DO ji = kideb , kiut 401 ztstemp(ji,jk) = t_s_b(ji,jk) 402 zeta_s(ji,jk) = rdt_ice / MAX(rhosn*cpic*zh_s(ji),epsi10) 403 403 END DO 404 404 END DO … … 443 443 DO numeq = nlay_s + 2, nlay_s + nlay_i 444 444 DO ji = kideb , kiut 445 layer= numeq - nlay_s - 1446 ztrid(ji,numeq,1) = - zeta_i(ji, layer)*zkappa_i(ji,layer-1)447 ztrid(ji,numeq,2) = 1.0 + zeta_i(ji, layer)*(zkappa_i(ji,layer-1) + &448 zkappa_i(ji, layer))449 ztrid(ji,numeq,3) = - zeta_i(ji, layer)*zkappa_i(ji,layer)450 zindterm(ji,numeq) = ztiold(ji, layer) + zeta_i(ji,layer)* &451 zradab_i(ji, layer)445 jk = numeq - nlay_s - 1 446 ztrid(ji,numeq,1) = - zeta_i(ji,jk)*zkappa_i(ji,jk-1) 447 ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,jk)*(zkappa_i(ji,jk-1) + & 448 zkappa_i(ji,jk)) 449 ztrid(ji,numeq,3) = - zeta_i(ji,jk)*zkappa_i(ji,jk) 450 zindterm(ji,numeq) = ztiold(ji,jk) + zeta_i(ji,jk)* & 451 zradab_i(ji,jk) 452 452 END DO 453 453 ENDDO … … 475 475 !!snow interior terms (bottom equation has the same form as the others) 476 476 DO numeq = 3, nlay_s + 1 477 layer= numeq - 1478 ztrid(ji,numeq,1) = - zeta_s(ji, layer)*zkappa_s(ji,layer-1)479 ztrid(ji,numeq,2) = 1.0 + zeta_s(ji, layer)*( zkappa_s(ji,layer-1) + &480 zkappa_s(ji, layer) )481 ztrid(ji,numeq,3) = - zeta_s(ji, layer)*zkappa_s(ji,layer)482 zindterm(ji,numeq) = ztsold(ji, layer) + zeta_s(ji,layer)* &483 zradab_s(ji, layer)477 jk = numeq - 1 478 ztrid(ji,numeq,1) = - zeta_s(ji,jk)*zkappa_s(ji,jk-1) 479 ztrid(ji,numeq,2) = 1.0 + zeta_s(ji,jk)*( zkappa_s(ji,jk-1) + & 480 zkappa_s(ji,jk) ) 481 ztrid(ji,numeq,3) = - zeta_s(ji,jk)*zkappa_s(ji,jk) 482 zindterm(ji,numeq) = ztsold(ji,jk) + zeta_s(ji,jk)* & 483 zradab_s(ji,jk) 484 484 END DO 485 485 … … 630 630 END DO 631 631 632 DO layer= minnumeqmin+1, maxnumeqmax633 DO ji = kideb , kiut 634 numeq = min(max(numeqmin(ji)+1, layer),numeqmax(ji))632 DO jk = minnumeqmin+1, maxnumeqmax 633 DO ji = kideb , kiut 634 numeq = min(max(numeqmin(ji)+1,jk),numeqmax(ji)) 635 635 zdiagbis(ji,numeq) = ztrid(ji,numeq,2) - ztrid(ji,numeq,1)* & 636 636 ztrid(ji,numeq-1,3)/zdiagbis(ji,numeq-1) … … 647 647 DO numeq = nlay_i + nlay_s + 1, nlay_s + 2, -1 648 648 DO ji = kideb , kiut 649 layer= numeq - nlay_s - 1650 t_i_b(ji, layer) = (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* &651 t_i_b(ji, layer+1))/zdiagbis(ji,numeq)649 jk = numeq - nlay_s - 1 650 t_i_b(ji,jk) = (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* & 651 t_i_b(ji,jk+1))/zdiagbis(ji,numeq) 652 652 END DO 653 653 END DO … … 679 679 END DO 680 680 681 DO layer= 1, nlay_s682 DO ji = kideb , kiut 683 t_s_b(ji, layer) = MAX( MIN( t_s_b(ji,layer), rtt ), 190._wp )684 zerrit(ji) = MAX(zerrit(ji),ABS(t_s_b(ji, layer) - ztstemp(ji,layer)))685 END DO 686 END DO 687 688 DO layer= 1, nlay_i689 DO ji = kideb , kiut 690 ztmelt_i = -tmut * s_i_b(ji, layer) + rtt691 t_i_b(ji, layer) = MAX(MIN(t_i_b(ji,layer),ztmelt_i), 190._wp)692 zerrit(ji) = MAX(zerrit(ji),ABS(t_i_b(ji, layer) - ztitemp(ji,layer)))681 DO jk = 1, nlay_s 682 DO ji = kideb , kiut 683 t_s_b(ji,jk) = MAX( MIN( t_s_b(ji,jk), rtt ), 190._wp ) 684 zerrit(ji) = MAX(zerrit(ji),ABS(t_s_b(ji,jk) - ztstemp(ji,jk))) 685 END DO 686 END DO 687 688 DO jk = 1, nlay_i 689 DO ji = kideb , kiut 690 ztmelt_i = -tmut * s_i_b(ji,jk) + rtt 691 t_i_b(ji,jk) = MAX(MIN(t_i_b(ji,jk),ztmelt_i), 190._wp) 692 zerrit(ji) = MAX(zerrit(ji),ABS(t_i_b(ji,jk) - ztitemp(ji,jk))) 693 693 END DO 694 694 END DO
Note: See TracChangeset
for help on using the changeset viewer.