Changeset 921 for trunk/NEMO/LIM_SRC_3/limthd_ent.F90
- Timestamp:
- 2008-05-13T10:28:52+02:00 (16 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/LIM_SRC_3/limthd_ent.F90 (modified) (36 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/LIM_SRC_3/limthd_ent.F90
r869 r921 46 46 CONTAINS 47 47 48 48 SUBROUTINE lim_thd_ent(kideb,kiut,jl) 49 49 !!------------------------------------------------------------------- 50 50 !! *** ROUTINE lim_thd_ent *** … … 135 135 REAL(wp), DIMENSION(jpij,0:jkmax+3) :: & 136 136 zhl0 ! old and new layer thicknesses 137 137 138 138 REAL(wp), DIMENSION(0:jkmax+3,0:jkmax+3) :: & 139 139 zrl01 … … 144 144 zqti_fin, zqts_fin 145 145 146 !------------------------------------------------------------------------------|146 !------------------------------------------------------------------------------| 147 147 148 148 zeps = 1.0d-20 … … 156 156 z_s(:,:) = 0.0 157 157 158 !159 !------------------------------------------------------------------------------|160 ! 1) Grid |161 !------------------------------------------------------------------------------|162 !158 ! 159 !------------------------------------------------------------------------------| 160 ! 1) Grid | 161 !------------------------------------------------------------------------------| 162 ! 163 163 nlays0 = nlay_s 164 164 nlayi0 = nlay_i … … 169 169 ENDDO 170 170 171 !172 !------------------------------------------------------------------------------|173 ! 2) Switches |174 !------------------------------------------------------------------------------|175 !171 ! 172 !------------------------------------------------------------------------------| 173 ! 2) Switches | 174 !------------------------------------------------------------------------------| 175 ! 176 176 ! 2.1 snind(ji), snswi(ji) 177 177 ! snow surface behaviour : computation of snind(ji)-snswi(ji) … … 181 181 ! 2 if 2nd layer is melting ... 182 182 DO ji = kideb, kiut 183 snind(ji) = 0184 zdeltah(ji) = 0.0183 snind(ji) = 0 184 zdeltah(ji) = 0.0 185 185 ENDDO !ji 186 186 187 187 DO jk = 1, nlays0 188 DO ji = kideb, kiut189 snind(ji) = jk * INT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)-zeps))) &190 191 zdeltah(ji)= zdeltah(ji) + zh_s(ji)192 END DO ! ji188 DO ji = kideb, kiut 189 snind(ji) = jk * INT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)-zeps))) & 190 + snind(ji) * (1 - INT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)-zeps)))) 191 zdeltah(ji)= zdeltah(ji) + zh_s(ji) 192 END DO ! ji 193 193 ENDDO ! jk 194 194 … … 198 198 snswi(ji) = MAX(0,INT(-dh_s_tot(ji)/MAX(zeps,ABS(dh_s_tot(ji))))) 199 199 ENDDO ! ji 200 200 201 201 ! 2.2 icsuind(ji), icsuswi(ji) 202 202 ! ice surface behaviour : computation of icsuind(ji)-icsuswi(ji) … … 206 206 ! 2 if 2nd layer is reached by melt ... 207 207 DO ji = kideb, kiut 208 icsuind(ji) = 0209 zdeltah(ji) = 0.0208 icsuind(ji) = 0 209 zdeltah(ji) = 0.0 210 210 ENDDO !ji 211 211 DO jk = 1, nlayi0 212 DO ji = kideb, kiut213 icsuind(ji) = jk * INT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)-zeps))) &214 215 zdeltah(ji) = zdeltah(ji) + zh_i(ji)216 END DO ! ji212 DO ji = kideb, kiut 213 icsuind(ji) = jk * INT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)-zeps))) & 214 + icsuind(ji) * (1 - INT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)-zeps)))) 215 zdeltah(ji) = zdeltah(ji) + zh_i(ji) 216 END DO ! ji 217 217 ENDDO !jk 218 218 … … 232 232 ! N+1 if all layers melt and that snow transforms into ice 233 233 DO ji = kideb, kiut 234 icboind(ji) = 0235 zdeltah(ji) = 0.0234 icboind(ji) = 0 235 zdeltah(ji) = 0.0 236 236 ENDDO 237 237 DO jk = nlayi0, 1, -1 238 DO ji = kideb, kiut239 icboind(ji) = (nlayi0+1-jk) &240 241 242 243 zdeltah(ji) = zdeltah(ji) + zh_i(ji)244 END DO238 DO ji = kideb, kiut 239 icboind(ji) = (nlayi0+1-jk) & 240 * INT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)-zeps))) & 241 + icboind(ji) & 242 * (1 - INT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)-zeps)))) 243 zdeltah(ji) = zdeltah(ji) + zh_i(ji) 244 END DO 245 245 ENDDO 246 246 … … 248 248 ! case of total ablation with remaining snow 249 249 IF ( ( ht_i_b(ji) .GT. zeps ) .AND. & 250 250 ( ht_i_b(ji) - dh_snowice(ji) .LT. zeps ) ) icboind(ji) = nlay_i + 1 251 251 END DO 252 252 … … 265 265 ! 2 if penultiem layer ... 266 266 DO ji = kideb, kiut 267 snicind(ji) = 0268 zdeltah(ji) = 0.0267 snicind(ji) = 0 268 zdeltah(ji) = 0.0 269 269 ENDDO 270 270 DO jk = nlays0, 1, -1 271 DO ji = kideb, kiut272 snicind(ji) = (nlays0+1-jk) &273 274 275 276 zdeltah(ji) = zdeltah(ji) + zh_s(ji)277 END DO271 DO ji = kideb, kiut 272 snicind(ji) = (nlays0+1-jk) & 273 * INT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji)-zeps))) & 274 + snicind(ji) & 275 * (1 - INT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji)-zeps)))) 276 zdeltah(ji) = zdeltah(ji) + zh_s(ji) 277 END DO 278 278 ENDDO 279 279 … … 282 282 ! 0 if not 283 283 DO ji = kideb, kiut 284 snicswi(ji) = MAX(0,INT(dh_snowice(ji)/MAX(zeps,ABS(dh_snowice(ji)))))285 ENDDO 286 287 !288 !------------------------------------------------------------------------------|289 ! 3) Snow redistribution |290 !------------------------------------------------------------------------------|291 !284 snicswi(ji) = MAX(0,INT(dh_snowice(ji)/MAX(zeps,ABS(dh_snowice(ji))))) 285 ENDDO 286 287 ! 288 !------------------------------------------------------------------------------| 289 ! 3) Snow redistribution | 290 !------------------------------------------------------------------------------| 291 ! 292 292 !------------- 293 293 ! Old profile … … 303 303 304 304 DO ji = kideb, kiut 305 nbot0(ji) = nlays0 + 1 - snind(ji) + ( 1. - snicind(ji) ) * &306 307 ! cotes of the top of the layers308 zm0(ji,0) = 0.0309 maxnbot0 = MAX ( maxnbot0 , nbot0(ji) )310 ENDDO 305 nbot0(ji) = nlays0 + 1 - snind(ji) + ( 1. - snicind(ji) ) * & 306 snicswi(ji) 307 ! cotes of the top of the layers 308 zm0(ji,0) = 0.0 309 maxnbot0 = MAX ( maxnbot0 , nbot0(ji) ) 310 ENDDO 311 311 IF( lk_mpp ) CALL mpp_max( maxnbot0, kcom=ncomm_ice ) 312 312 313 313 DO jk = 1, maxnbot0 314 DO ji = kideb, kiut315 !change316 limsum = ( 1 - snswi(ji) ) * ( jk - 1 ) + &317 318 limsum = MIN( limsum , nlay_s )319 zm0(ji,jk) = dh_s_tot(ji) + zh_s(ji) * limsum320 END DO314 DO ji = kideb, kiut 315 !change 316 limsum = ( 1 - snswi(ji) ) * ( jk - 1 ) + & 317 snswi(ji) * ( jk + snind(ji) - 1 ) 318 limsum = MIN( limsum , nlay_s ) 319 zm0(ji,jk) = dh_s_tot(ji) + zh_s(ji) * limsum 320 END DO 321 321 ENDDO 322 322 323 323 DO ji = kideb, kiut 324 324 zm0(ji,nbot0(ji)) = dh_s_tot(ji) - snicswi(ji) * dh_snowice(ji) + & 325 325 zh_s(ji) * nlays0 326 326 zm0(ji,1) = dh_s_tot(ji) * (1 -snswi(ji) ) + & 327 327 snswi(ji) * zm0(ji,1) 328 328 ENDDO 329 329 330 330 DO jk = ntop0, maxnbot0 331 DO ji = kideb, kiut332 ! layer thickness333 zthick0(ji,jk) = zm0(ji,jk) - zm0(ji,jk-1)334 END DO331 DO ji = kideb, kiut 332 ! layer thickness 333 zthick0(ji,jk) = zm0(ji,jk) - zm0(ji,jk-1) 334 END DO 335 335 ENDDO 336 336 337 337 zqts_in(:) = 0.0 338 339 DO ji = kideb, kiut 340 ! layer heat content341 qm0(ji,1) = rhosn * ( cpic * ( rtt - ( 1. - snswi(ji) ) * ( tatm_ice_1d(ji) ) &342 343 344 zqts_in(ji) = zqts_in(ji) + qm0(ji,1)338 339 DO ji = kideb, kiut 340 ! layer heat content 341 qm0(ji,1) = rhosn * ( cpic * ( rtt - ( 1. - snswi(ji) ) * ( tatm_ice_1d(ji) ) & 342 - snswi(ji) * t_s_b(ji,1) ) & 343 + lfus ) * zthick0(ji,1) 344 zqts_in(ji) = zqts_in(ji) + qm0(ji,1) 345 345 ENDDO 346 346 347 347 DO jk = 2, maxnbot0 348 DO ji = kideb, kiut349 limsum = ( 1 - snswi(ji) ) * ( jk - 1 ) + &350 351 limsum = MIN( limsum , nlay_s )352 qm0(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,limsum) ) + lfus ) &353 354 zswitch = 1.0 - MAX (0.0, SIGN ( 1.0, zeps - ht_s_b(ji) ) )355 zqts_in(ji) = zqts_in(ji) + ( 1. - snswi(ji) ) * qm0(ji,jk) * zswitch356 END DO ! jk348 DO ji = kideb, kiut 349 limsum = ( 1 - snswi(ji) ) * ( jk - 1 ) + & 350 snswi(ji) * ( jk + snind(ji) - 1 ) 351 limsum = MIN( limsum , nlay_s ) 352 qm0(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,limsum) ) + lfus ) & 353 * zthick0(ji,jk) 354 zswitch = 1.0 - MAX (0.0, SIGN ( 1.0, zeps - ht_s_b(ji) ) ) 355 zqts_in(ji) = zqts_in(ji) + ( 1. - snswi(ji) ) * qm0(ji,jk) * zswitch 356 END DO ! jk 357 357 ENDDO ! ji 358 358 … … 362 362 ! zqsnow, enthalpy of the flooded snow 363 363 DO ji = kideb, kiut 364 zqsnow(ji) = rhosn*lfus365 zdeltah(ji) = 0.0364 zqsnow(ji) = rhosn*lfus 365 zdeltah(ji) = 0.0 366 366 ENDDO 367 367 368 368 DO jk = nlays0, 1, -1 369 DO ji = kideb, kiut370 zhsnow = MAX(0.0,dh_snowice(ji)-zdeltah(ji))371 zqsnow(ji) = zqsnow(ji) + &372 373 zdeltah(ji) = zdeltah(ji) + zh_s(ji)374 END DO369 DO ji = kideb, kiut 370 zhsnow = MAX(0.0,dh_snowice(ji)-zdeltah(ji)) 371 zqsnow(ji) = zqsnow(ji) + & 372 rhosn*cpic*(rtt-t_s_b(ji,jk)) 373 zdeltah(ji) = zdeltah(ji) + zh_s(ji) 374 END DO 375 375 ENDDO 376 376 … … 398 398 399 399 DO jk = 1, nlay_s 400 DO ji = kideb, kiut401 z_s(ji,jk) = zh_s(ji) * jk402 END DO400 DO ji = kideb, kiut 401 z_s(ji,jk) = zh_s(ji) * jk 402 END DO 403 403 ENDDO 404 404 … … 407 407 !----------------- 408 408 DO layer0 = ntop0, maxnbot0 409 DO ji = kideb, kiut410 zhl0(ji,layer0) = zm0(ji,layer0) - zm0(ji,layer0-1)411 END DO409 DO ji = kideb, kiut 410 zhl0(ji,layer0) = zm0(ji,layer0) - zm0(ji,layer0-1) 411 END DO 412 412 ENDDO 413 413 414 414 DO layer1 = ntop1, nbot1 415 DO ji = kideb, kiut416 q_s_b(ji,layer1)= 0.0417 END DO415 DO ji = kideb, kiut 416 q_s_b(ji,layer1)= 0.0 417 END DO 418 418 ENDDO 419 419 … … 422 422 !---------------- 423 423 DO layer0 = ntop0, maxnbot0 424 DO layer1 = ntop1, nbot1425 DO ji = kideb, kiut426 zrl01(layer1,layer0) = MAX(0.0,( MIN(zm0(ji,layer0),z_s(ji,layer1)) &427 - MAX(zm0(ji,layer0-1), z_s(ji,layer1-1)))/MAX(zhl0(ji,layer0),epsi10))428 q_s_b(ji,layer1) = q_s_b(ji,layer1) + zrl01(layer1,layer0)*qm0(ji,layer0) &429 430 END DO431 END DO424 DO layer1 = ntop1, nbot1 425 DO ji = kideb, kiut 426 zrl01(layer1,layer0) = MAX(0.0,( MIN(zm0(ji,layer0),z_s(ji,layer1)) & 427 - MAX(zm0(ji,layer0-1), z_s(ji,layer1-1)))/MAX(zhl0(ji,layer0),epsi10)) 428 q_s_b(ji,layer1) = q_s_b(ji,layer1) + zrl01(layer1,layer0)*qm0(ji,layer0) & 429 * MAX(0.0,SIGN(1.0,nbot0(ji)-layer0+zeps)) 430 END DO 431 END DO 432 432 ENDDO 433 433 … … 441 441 442 442 IF ( con_i ) THEN 443 DO ji = kideb, kiut444 IF ( ABS ( zqts_in(ji) - zqts_fin(ji) ) / rdt_ice .GT. 1.0e-6 ) THEN445 zji = MOD( npb(ji) - 1, jpi ) + 1446 zjj = ( npb(ji) - 1 ) / jpi + 1447 WRITE(numout,*) ' violation of heat conservation : ', &448 449 WRITE(numout,*) ' ji, jj : ', zji, zjj450 WRITE(numout,*) ' ht_s_b : ', ht_s_b(ji)451 WRITE(numout,*) ' zqts_in : ', zqts_in(ji) / rdt_ice452 WRITE(numout,*) ' zqts_fin : ', zqts_fin(ji) / rdt_ice453 WRITE(numout,*) ' dh_snowice : ', dh_snowice(ji)454 WRITE(numout,*) ' dh_s_tot : ', dh_s_tot(ji)455 WRITE(numout,*) ' snswi : ', snswi(ji)456 ENDIF457 END DO443 DO ji = kideb, kiut 444 IF ( ABS ( zqts_in(ji) - zqts_fin(ji) ) / rdt_ice .GT. 1.0e-6 ) THEN 445 zji = MOD( npb(ji) - 1, jpi ) + 1 446 zjj = ( npb(ji) - 1 ) / jpi + 1 447 WRITE(numout,*) ' violation of heat conservation : ', & 448 ABS ( zqts_in(ji) - zqts_fin(ji) ) / rdt_ice 449 WRITE(numout,*) ' ji, jj : ', zji, zjj 450 WRITE(numout,*) ' ht_s_b : ', ht_s_b(ji) 451 WRITE(numout,*) ' zqts_in : ', zqts_in(ji) / rdt_ice 452 WRITE(numout,*) ' zqts_fin : ', zqts_fin(ji) / rdt_ice 453 WRITE(numout,*) ' dh_snowice : ', dh_snowice(ji) 454 WRITE(numout,*) ' dh_s_tot : ', dh_s_tot(ji) 455 WRITE(numout,*) ' snswi : ', snswi(ji) 456 ENDIF 457 END DO 458 458 ENDIF 459 459 … … 473 473 zfac2 = lfus / cpic 474 474 DO jk = 1, nlay_s 475 DO ji = kideb, kiut476 zswitch = MAX ( 0.0 , SIGN ( 1.0, zeps - ht_s_b(ji) ) )477 t_s_b(ji,jk) = rtt &478 479 480 END DO481 ENDDO 482 !483 !------------------------------------------------------------------------------|484 ! 4) Ice redistribution |485 !------------------------------------------------------------------------------|486 !475 DO ji = kideb, kiut 476 zswitch = MAX ( 0.0 , SIGN ( 1.0, zeps - ht_s_b(ji) ) ) 477 t_s_b(ji,jk) = rtt & 478 + ( 1.0 - zswitch ) * & 479 ( - zfac1 * q_s_b(ji,jk) + zfac2 ) 480 END DO 481 ENDDO 482 ! 483 !------------------------------------------------------------------------------| 484 ! 4) Ice redistribution | 485 !------------------------------------------------------------------------------| 486 ! 487 487 !------------- 488 488 ! OLD PROFILE … … 496 496 497 497 DO ji = kideb, kiut 498 ! reference number of the bottommost layer498 ! reference number of the bottommost layer 499 499 nbot0(ji) = MAX( 1 , MIN( nlayi0 + ( 1 - icboind(ji) ) + & 500 501 500 ( 1 - icsuind(ji) ) * icsuswi(ji) + snicswi(ji) , & 501 nlay_i + 2 ) ) 502 502 ! maximum reference number of the bottommost layer over all domain 503 503 maxnbot0 = MAX( maxnbot0 , nbot0(ji) ) … … 508 508 !------------------------- 509 509 zm0(:,0) = 0.0 510 510 511 511 DO jk = 1, maxnbot0 512 512 DO ji = kideb, kiut … … 515 515 ! limsum is the real ice layer number corresponding to present jk 516 516 limsum = ( (icsuswi(ji)*(icsuind(ji)+jk-1) + & 517 517 (1-icsuswi(ji))*jk))*(1-snicswi(ji)) + (jk-1)*snicswi(ji) 518 518 zm0(ji,jk)= icsuswi(ji)*dh_i_surf(ji) + snicswi(ji)*dh_snowice(ji) & 519 519 + limsum * zh_i(ji) 520 520 END DO 521 521 ENDDO … … 523 523 DO ji = kideb, kiut 524 524 zm0(ji,nbot0(ji)) = icsuswi(ji)*dh_i_surf(ji) + snicswi(ji)*dh_snowice(ji) + dh_i_bott(ji) & 525 525 + zh_i(ji) * nlayi0 526 526 zm0(ji,1) = snicswi(ji)*dh_snowice(ji) + (1-snicswi(ji))*zm0(ji,1) 527 527 ENDDO … … 531 531 !----------------------------- 532 532 DO jk = ntop0, maxnbot0 533 DO ji = kideb, kiut534 zthick0(ji,jk) = zm0(ji,jk) - zm0(ji,jk-1)535 END DO533 DO ji = kideb, kiut 534 zthick0(ji,jk) = zm0(ji,jk) - zm0(ji,jk-1) 535 END DO 536 536 ENDDO 537 537 … … 545 545 DO ji = kideb, kiut 546 546 limsum = MAX(1,MIN(snicswi(ji)*(jk-1) + icsuswi(ji)*(jk-1+icsuind(ji)) + & 547 547 (1-icsuswi(ji))*(1-snicswi(ji))*jk,nlay_i)) 548 548 ztmelts = -tmut * s_i_b(ji,limsum) + rtt 549 549 qm0(ji,jk) = rhoic * ( cpic * (ztmelts-t_i_b(ji,limsum)) + lfus * ( 1.0-(ztmelts-rtt)/ & 550 551 550 MIN((t_i_b(ji,limsum)-rtt),-zeps) ) - rcp*(ztmelts-rtt) ) & 551 * zthick0(ji,jk) 552 552 END DO 553 553 ENDDO … … 557 557 !---------------------------- 558 558 DO ji = kideb, kiut 559 ztmelts = ( 1.0 - icboswi(ji) ) * (-tmut * s_i_b(ji,nlayi0)) & ! case of melting ice560 561 562 563 ! bottom formation temperature564 ztform = t_i_b(ji,nlay_i)565 IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) ztform = t_bo_b(ji)566 qm0(ji,nbot0(ji)) = ( 1.0 - icboswi(ji) )*qm0(ji,nbot0(ji)) & ! case of melting ice567 568 569 570 571 572 559 ztmelts = ( 1.0 - icboswi(ji) ) * (-tmut * s_i_b(ji,nlayi0)) & ! case of melting ice 560 + icboswi(ji) * (-tmut * s_i_new(ji)) & ! case of forming ice 561 + rtt ! this temperature is in Celsius 562 563 ! bottom formation temperature 564 ztform = t_i_b(ji,nlay_i) 565 IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) ztform = t_bo_b(ji) 566 qm0(ji,nbot0(ji)) = ( 1.0 - icboswi(ji) )*qm0(ji,nbot0(ji)) & ! case of melting ice 567 + icboswi(ji) * & ! case of forming ice 568 rhoic*( cpic*(ztmelts-ztform) & 569 + lfus *( 1.0-(ztmelts-rtt)/ & 570 MIN ( (ztform-rtt) , - epsi10 ) ) & 571 - rcp*(ztmelts-rtt) ) & 572 *zthick0(ji,nbot0(ji)) 573 573 ENDDO 574 574 … … 579 579 ! energy of the flooding seawater 580 580 zqsnic = rau0 * rcp * ( rtt - t_bo_b(ji) ) * dh_snowice(ji) * & 581 581 (rhoic - rhosn) / rhoic * snicswi(ji) ! generally positive 582 582 ! Heat conservation diagnostic 583 583 qt_i_in(ji,jl) = qt_i_in(ji,jl) + zqsnic … … 593 593 594 594 DO jk = ntop0, maxnbot0 595 DO ji = kideb, kiut596 ! Heat conservation597 zqti_in(ji) = zqti_in(ji) + qm0(ji,jk) &598 599 600 END DO595 DO ji = kideb, kiut 596 ! Heat conservation 597 zqti_in(ji) = zqti_in(ji) + qm0(ji,jk) & 598 * MAX( 0.0 , SIGN(1.0,ht_i_b(ji)-zeps6+zeps) ) & 599 * MAX( 0.0 , SIGN( 1. , nbot0(ji) - jk + zeps ) ) 600 END DO 601 601 ENDDO 602 602 … … 616 616 !------------------ 617 617 DO ji = kideb, kiut 618 zh_i(ji) = ht_i_b(ji) / nlay_i618 zh_i(ji) = ht_i_b(ji) / nlay_i 619 619 ENDDO 620 620 … … 624 624 z_i(:,0) = 0.0 625 625 DO jk = 1, nlay_i 626 DO ji = kideb, kiut627 z_i(ji,jk) = zh_i(ji) * jk628 END DO626 DO ji = kideb, kiut 627 z_i(ji,jk) = zh_i(ji) * jk 628 END DO 629 629 ENDDO 630 630 631 631 !--thicknesses of the layers 632 632 DO layer0 = ntop0, maxnbot0 633 DO ji = kideb, kiut634 zhl0(ji,layer0) = zm0(ji,layer0) - zm0(ji,layer0-1) !thicknesses of the layers635 END DO633 DO ji = kideb, kiut 634 zhl0(ji,layer0) = zm0(ji,layer0) - zm0(ji,layer0-1) !thicknesses of the layers 635 END DO 636 636 ENDDO 637 637 … … 642 642 q_i_b(:,:) = 0.0 643 643 DO layer0 = ntop0, maxnbot0 644 DO layer1 = ntop1, nbot1645 DO ji = kideb, kiut646 zrl01(layer1,layer0) = MAX(0.0,( MIN(zm0(ji,layer0),z_i(ji,layer1)) &647 - MAX(zm0(ji,layer0-1), z_i(ji,layer1-1)))/MAX(zhl0(ji,layer0),epsi10))648 q_i_b(ji,layer1) = q_i_b(ji,layer1) &649 650 651 652 END DO653 END DO644 DO layer1 = ntop1, nbot1 645 DO ji = kideb, kiut 646 zrl01(layer1,layer0) = MAX(0.0,( MIN(zm0(ji,layer0),z_i(ji,layer1)) & 647 - MAX(zm0(ji,layer0-1), z_i(ji,layer1-1)))/MAX(zhl0(ji,layer0),epsi10)) 648 q_i_b(ji,layer1) = q_i_b(ji,layer1) & 649 + zrl01(layer1,layer0)*qm0(ji,layer0) & 650 * MAX(0.0,SIGN(1.0,ht_i_b(ji)-zeps6+zeps)) & 651 * MAX(0.0,SIGN(1.0,nbot0(ji)-layer0+zeps)) 652 END DO 653 END DO 654 654 ENDDO 655 655 … … 663 663 END DO 664 664 END DO 665 !665 ! 666 666 DO ji = kideb, kiut 667 667 IF ( ABS ( zqti_in(ji) - zqti_fin(ji) ) / rdt_ice .GT. 1.0e-6 ) THEN … … 669 669 zjj = ( npb(ji) - 1 ) / jpi + 1 670 670 WRITE(numout,*) ' violation of heat conservation : ', & 671 671 ABS ( zqti_in(ji) - zqti_fin(ji) ) / rdt_ice 672 672 WRITE(numout,*) ' ji, jj : ', zji, zjj 673 673 WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji) … … 700 700 END DO 701 701 702 !703 !------------------------------------------------------------------------------|704 ! 5) Update salinity and recover temperature |705 !------------------------------------------------------------------------------|706 !702 ! 703 !------------------------------------------------------------------------------| 704 ! 5) Update salinity and recover temperature | 705 !------------------------------------------------------------------------------| 706 ! 707 707 ! Update salinity (basal entrapment, snow ice formation) 708 708 DO ji = kideb, kiut 709 709 sm_i_b(ji) = sm_i_b(ji) & 710 710 + dsm_i_se_1d(ji) + dsm_i_si_1d(ji) 711 711 END DO !ji 712 712 … … 720 720 zaaa = cpic 721 721 zbbb = ( rcp - cpic ) * ( ztmelts - rtt ) + & 722 722 q_i_b(ji,jk) / rhoic - lfus 723 723 zccc = lfus * ( ztmelts - rtt ) 724 724 zdiscrim = SQRT( MAX(zbbb*zbbb - 4.0*zaaa*zccc,0.0) ) 725 725 t_i_b(ji,jk) = rtt - ( zbbb + zdiscrim ) / & 726 726 ( 2.0 *zaaa ) 727 727 END DO !ji 728 728 729 729 END DO !jk 730 730 731 731 END SUBROUTINE lim_thd_ent 732 732 733 733 #else … … 740 740 END SUBROUTINE lim_thd_ent 741 741 #endif 742 742 END MODULE limthd_ent
Note: See TracChangeset
for help on using the changeset viewer.
