- Timestamp:
- 2013-11-07T11:01:27+01:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90
r3625 r4161 44 44 45 45 !!---------------------------------------------------------------------- 46 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)46 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 47 47 !! $Id$ 48 48 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 75 75 76 76 INTEGER :: ji,jk ! dummy loop indices 77 INTEGER :: zji, zjj , & ! dummy indices77 INTEGER :: ii, ij , & ! dummy indices 78 78 ntop0 , & ! old layer top index 79 79 nbot1 , & ! new layer bottom index … … 145 145 146 146 DO ji = kideb, kiut 147 zh_i(ji) = old_ht_i_b(ji) / nlay_i148 zh_s(ji) = old_ht_s_b(ji) / nlay_s147 zh_i(ji) = old_ht_i_b(ji) / REAL( nlay_i ) 148 zh_s(ji) = old_ht_s_b(ji) / REAL( nlay_s ) 149 149 END DO 150 150 … … 166 166 DO jk = 1, nlays0 167 167 DO ji = kideb, kiut 168 snind(ji) = jk * INT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)-epsi20))) &169 + snind(ji) * (1 - INT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)-epsi20))))168 snind(ji) = jk * NINT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)))) & 169 + snind(ji) * (1 - NINT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji))))) 170 170 zdeltah(ji)= zdeltah(ji) + zh_s(ji) 171 171 END DO ! ji … … 175 175 ! 0 if not 176 176 DO ji = kideb, kiut 177 snswi(ji) = MAX(0, INT(-dh_s_tot(ji)/MAX(epsi20,ABS(dh_s_tot(ji)))))177 snswi(ji) = MAX(0,NINT(-dh_s_tot(ji)/MAX(epsi20,ABS(dh_s_tot(ji))))) 178 178 END DO ! ji 179 179 … … 190 190 DO jk = 1, nlayi0 191 191 DO ji = kideb, kiut 192 icsuind(ji) = jk * INT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)-epsi20))) &193 + icsuind(ji) * (1 - INT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)-epsi20))))192 icsuind(ji) = jk * NINT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)))) & 193 + icsuind(ji) * (1 - NINT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji))))) 194 194 zdeltah(ji) = zdeltah(ji) + zh_i(ji) 195 195 END DO ! ji … … 200 200 ! 0 if not 201 201 DO ji = kideb, kiut 202 icsuswi(ji) = MAX(0, INT(-dh_i_surf(ji)/MAX(epsi20 , ABS(dh_i_surf(ji)) ) ) )202 icsuswi(ji) = MAX(0,NINT(-dh_i_surf(ji)/MAX(epsi20 , ABS(dh_i_surf(ji)) ) ) ) 203 203 ENDDO 204 204 … … 216 216 DO jk = nlayi0, 1, -1 217 217 DO ji = kideb, kiut 218 icboind(ji) = (nlayi0+1-jk) * INT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)-epsi20))) &219 & + icboind(ji) * (1 - INT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)-epsi20))))218 icboind(ji) = (nlayi0+1-jk) * NINT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)))) & 219 & + icboind(ji) * (1 - NINT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji))))) 220 220 zdeltah(ji) = zdeltah(ji) + zh_i(ji) 221 221 END DO … … 232 232 ! 0 if ablation is on the way 233 233 DO ji = kideb, kiut 234 icboswi(ji) = MAX(0, INT(dh_i_bott(ji) / MAX(epsi20,ABS(dh_i_bott(ji)))))234 icboswi(ji) = MAX(0,NINT(dh_i_bott(ji) / MAX(epsi20,ABS(dh_i_bott(ji))))) 235 235 END DO 236 236 … … 248 248 DO ji = kideb, kiut 249 249 snicind(ji) = (nlays0+1-jk) & 250 * INT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji)-epsi20))) + snicind(ji) &251 * (1 - INT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji)-epsi20))))250 * NINT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji)))) + snicind(ji) & 251 * (1 - NINT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji))))) 252 252 zdeltah(ji) = zdeltah(ji) + zh_s(ji) 253 253 END DO … … 258 258 ! 0 if not 259 259 DO ji = kideb, kiut 260 snicswi(ji) = MAX(0, INT(dh_snowice(ji)/MAX(epsi20,ABS(dh_snowice(ji)))))260 snicswi(ji) = MAX(0,NINT(dh_snowice(ji)/MAX(epsi20,ABS(dh_snowice(ji))))) 261 261 ENDDO 262 262 … … 279 279 280 280 DO ji = kideb, kiut 281 nbot0(ji) = nlays0 + 1 - snind(ji) + ( 1 .- snicind(ji) ) * snicswi(ji)281 nbot0(ji) = nlays0 + 1 - snind(ji) + ( 1 - snicind(ji) ) * snicswi(ji) 282 282 ! cotes of the top of the layers 283 283 zm0(ji,0) = 0._wp … … 291 291 limsum = ( 1 - snswi(ji) ) * ( jk - 1 ) + snswi(ji) * ( jk + snind(ji) - 1 ) 292 292 limsum = MIN( limsum , nlay_s ) 293 zm0(ji,jk) = dh_s_tot(ji) + zh_s(ji) * limsum294 END DO 295 END DO 296 297 DO ji = kideb, kiut 298 zm0(ji,nbot0(ji)) = dh_s_tot(ji) - snicswi(ji) * dh_snowice(ji) + zh_s(ji) * nlays0299 zm0(ji,1) = dh_s_tot(ji) * (1 -snswi(ji) ) + snswi(ji) * zm0(ji,1)293 zm0(ji,jk) = dh_s_tot(ji) + zh_s(ji) * REAL( limsum ) 294 END DO 295 END DO 296 297 DO ji = kideb, kiut 298 zm0(ji,nbot0(ji)) = dh_s_tot(ji) - REAL( snicswi(ji) ) * dh_snowice(ji) + zh_s(ji) * REAL( nlays0 ) 299 zm0(ji,1) = dh_s_tot(ji) * REAL( 1 - snswi(ji) ) + REAL( snswi(ji) ) * zm0(ji,1) 300 300 END DO 301 301 … … 309 309 310 310 DO ji = kideb, kiut ! layer heat content 311 qm0 (ji,1) = rhosn * ( cpic * ( rtt - ( 1.- snswi(ji) ) * tatm_ice_1d(ji) &312 & - snswi(ji)* t_s_b (ji,1) ) &311 qm0 (ji,1) = rhosn * ( cpic * ( rtt - REAL( 1 - snswi(ji) ) * tatm_ice_1d(ji) & 312 & - REAL( snswi(ji) ) * t_s_b (ji,1) ) & 313 313 & + lfus ) * zthick0(ji,1) 314 314 zqts_in(ji) = zqts_in(ji) + qm0(ji,1) … … 320 320 limsum = MIN( limsum , nlay_s ) 321 321 qm0(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,limsum) ) + lfus ) * zthick0(ji,jk) 322 zswitch = 1.0 - MAX (0.0, SIGN ( 1.0, epsi20- ht_s_b(ji) ) )323 zqts_in(ji) = zqts_in(ji) + ( 1.- snswi(ji) ) * qm0(ji,jk) * zswitch322 zswitch = 1.0 - MAX (0.0, SIGN ( 1.0, - ht_s_b(ji) ) ) 323 zqts_in(ji) = zqts_in(ji) + REAL( 1 - snswi(ji) ) * qm0(ji,jk) * zswitch 324 324 END DO ! jk 325 325 END DO ! ji … … 360 360 !------------------- 361 361 DO ji = kideb, kiut 362 zh_s(ji) = ht_s_b(ji) / nlay_s362 zh_s(ji) = ht_s_b(ji) / REAL( nlay_s ) 363 363 z_s(ji,0) = 0._wp 364 364 ENDDO … … 366 366 DO jk = 1, nlay_s 367 367 DO ji = kideb, kiut 368 z_s(ji,jk) = zh_s(ji) * jk368 z_s(ji,jk) = zh_s(ji) * REAL( jk ) 369 369 END DO 370 370 END DO … … 394 394 & - MAX(zm0(ji,layer0-1), z_s(ji,layer1-1))) / MAX(zhl0(ji,layer0),epsi10)) 395 395 q_s_b(ji,layer1) = q_s_b(ji,layer1) + zrl01(layer1,layer0)*qm0(ji,layer0) & 396 & * MAX(0.0,SIGN(1.0, nbot0(ji)-layer0+epsi20))396 & * MAX(0.0,SIGN(1.0,REAL(nbot0(ji)-layer0))) 397 397 END DO 398 398 END DO … … 410 410 DO ji = kideb, kiut 411 411 IF ( ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice > 1.0e-6 ) THEN 412 zji = MOD( npb(ji) - 1, jpi ) + 1 413 zjj = ( npb(ji) - 1 ) / jpi + 1 414 WRITE(numout,*) ' violation of heat conservation : ', & 415 ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice 416 WRITE(numout,*) ' ji, jj : ', zji, zjj 412 ii = MOD( npb(ji) - 1, jpi ) + 1 413 ij = ( npb(ji) - 1 ) / jpi + 1 414 WRITE(numout,*) ' violation of heat conservation : ', ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice 415 WRITE(numout,*) ' ji, jj : ', ii, ij 417 416 WRITE(numout,*) ' ht_s_b : ', ht_s_b(ji) 418 417 WRITE(numout,*) ' zqts_in : ', zqts_in (ji) * r1_rdtice … … 441 440 DO jk = 1, nlay_s 442 441 DO ji = kideb, kiut 443 zswitch = MAX ( 0.0 , SIGN ( 1.0, epsi20- ht_s_b(ji) ) )442 zswitch = MAX ( 0.0 , SIGN ( 1.0, - ht_s_b(ji) ) ) 444 443 t_s_b(ji,jk) = rtt + ( 1.0 - zswitch ) * ( - zfac1 * q_s_b(ji,jk) + zfac2 ) 445 444 END DO … … 480 479 limsum = ( (icsuswi(ji)*(icsuind(ji)+jk-1) + & 481 480 (1-icsuswi(ji))*jk))*(1-snicswi(ji)) + (jk-1)*snicswi(ji) 482 zm0(ji,jk)= icsuswi(ji)*dh_i_surf(ji) + snicswi(ji)*dh_snowice(ji) &483 + limsum* zh_i(ji)484 END DO 485 END DO 486 487 DO ji = kideb, kiut 488 zm0(ji,nbot0(ji)) = icsuswi(ji)*dh_i_surf(ji) + snicswi(ji)*dh_snowice(ji) + dh_i_bott(ji) &489 + zh_i(ji) * nlayi0490 zm0(ji,1) = snicswi(ji)*dh_snowice(ji) +(1-snicswi(ji))*zm0(ji,1)481 zm0(ji,jk)= REAL(icsuswi(ji))*dh_i_surf(ji) + REAL(snicswi(ji))*dh_snowice(ji) & 482 + REAL(limsum) * zh_i(ji) 483 END DO 484 END DO 485 486 DO ji = kideb, kiut 487 zm0(ji,nbot0(ji)) = REAL(icsuswi(ji))*dh_i_surf(ji) + REAL(snicswi(ji))*dh_snowice(ji) + dh_i_bott(ji) & 488 + zh_i(ji) * REAL(nlayi0) 489 zm0(ji,1) = REAL(snicswi(ji))*dh_snowice(ji) + REAL(1-snicswi(ji))*zm0(ji,1) 491 490 END DO 492 491 … … 521 520 !---------------------------- 522 521 DO ji = kideb, kiut 523 ztmelts = ( 1.0- icboswi(ji) ) * (-tmut * s_i_b (ji,nlayi0) ) & ! case of melting ice524 & + icboswi(ji)* (-tmut * s_i_new(ji) ) & ! case of forming ice522 ztmelts = REAL( 1 - icboswi(ji) ) * (-tmut * s_i_b (ji,nlayi0) ) & ! case of melting ice 523 & + REAL( icboswi(ji) ) * (-tmut * s_i_new(ji) ) & ! case of forming ice 525 524 & + rtt ! in Kelvin 526 525 … … 528 527 ztform = t_i_b(ji,nlay_i) 529 528 IF( num_sal == 2 ) ztform = t_bo_b(ji) 530 qm0(ji,nbot0(ji)) = ( 1.0- icboswi(ji) )*qm0(ji,nbot0(ji)) & ! case of melting ice531 & + icboswi(ji) * rhoic * ( cpic*(ztmelts-ztform) & ! case of forming ice529 qm0(ji,nbot0(ji)) = REAL( 1 - icboswi(ji) )*qm0(ji,nbot0(ji)) & ! case of melting ice 530 & + REAL( icboswi(ji) ) * rhoic * ( cpic*(ztmelts-ztform) & ! case of forming ice 532 531 + lfus *( 1.0-(ztmelts-rtt) / MIN ( (ztform-rtt) , - epsi10 ) ) & 533 532 - rcp*(ztmelts-rtt) ) * zthick0(ji,nbot0(ji) ) … … 540 539 ! energy of the flooding seawater 541 540 zqsnic = rau0 * rcp * ( rtt - t_bo_b(ji) ) * dh_snowice(ji) * & 542 (rhoic - rhosn) / rhoic * snicswi(ji) ! generally positive541 (rhoic - rhosn) / rhoic * REAL(snicswi(ji)) ! generally positive 543 542 ! Heat conservation diagnostic 544 543 qt_i_in(ji,jl) = qt_i_in(ji,jl) + zqsnic … … 549 548 ! = enthalpy of snow + enthalpy of frozen water 550 549 zqsnic = zqsnow(ji) + zqsnic 551 qm0(ji,1) = snicswi(ji) * zqsnic +( 1 - snicswi(ji) ) * qm0(ji,1)550 qm0(ji,1) = REAL(snicswi(ji)) * zqsnic + REAL( 1 - snicswi(ji) ) * qm0(ji,1) 552 551 553 552 END DO ! ji … … 556 555 DO ji = kideb, kiut 557 556 ! Heat conservation 558 zqti_in(ji) = zqti_in(ji) + qm0(ji,jk) * MAX( 0.0 , SIGN(1.0,ht_i_b(ji)-epsi06 +epsi20) ) &559 & * MAX( 0.0 , SIGN( 1. , nbot0(ji) - jk + epsi20) )557 zqti_in(ji) = zqti_in(ji) + qm0(ji,jk) * MAX( 0.0 , SIGN(1.0,ht_i_b(ji)-epsi06) ) & 558 & * MAX( 0.0 , SIGN( 1. , REAL(nbot0(ji) - jk) ) ) 560 559 END DO 561 560 END DO … … 575 574 !------------------ 576 575 DO ji = kideb, kiut 577 zh_i(ji) = ht_i_b(ji) / nlay_i576 zh_i(ji) = ht_i_b(ji) / REAL( nlay_i ) 578 577 ENDDO 579 578 … … 606 605 q_i_b(ji,layer1) = q_i_b(ji,layer1) & 607 606 + zrl01(layer1,layer0)*qm0(ji,layer0) & 608 * MAX(0.0,SIGN(1.0,ht_i_b(ji)-epsi06 +epsi20)) &609 * MAX(0.0,SIGN(1.0, nbot0(ji)-layer0+epsi20))607 * MAX(0.0,SIGN(1.0,ht_i_b(ji)-epsi06)) & 608 * MAX(0.0,SIGN(1.0,REAL(nbot0(ji)-layer0))) 610 609 END DO 611 610 END DO … … 622 621 END DO 623 622 ! 624 DO ji = kideb, kiut 625 IF ( ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice > 1.0e-6 ) THEN 626 zji = MOD( npb(ji) - 1, jpi ) + 1 627 zjj = ( npb(ji) - 1 ) / jpi + 1 628 WRITE(numout,*) ' violation of heat conservation : ', ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice 629 WRITE(numout,*) ' ji, jj : ', zji, zjj 630 WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji) 631 WRITE(numout,*) ' zqti_in : ', zqti_in (ji) * r1_rdtice 632 WRITE(numout,*) ' zqti_fin : ', zqti_fin(ji) * r1_rdtice 633 WRITE(numout,*) ' dh_i_bott: ', dh_i_bott(ji) 634 WRITE(numout,*) ' dh_i_surf: ', dh_i_surf(ji) 635 WRITE(numout,*) ' dh_snowice:', dh_snowice(ji) 636 WRITE(numout,*) ' icsuswi : ', icsuswi(ji) 637 WRITE(numout,*) ' icboswi : ', icboswi(ji) 638 WRITE(numout,*) ' snicswi : ', snicswi(ji) 639 ENDIF 640 END DO 623 IF ( con_i ) THEN 624 DO ji = kideb, kiut 625 IF ( ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice > 1.0e-6 ) THEN 626 ii = MOD( npb(ji) - 1, jpi ) + 1 627 ij = ( npb(ji) - 1 ) / jpi + 1 628 WRITE(numout,*) ' violation of heat conservation : ', ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice 629 WRITE(numout,*) ' ji, jj : ', ii, ij 630 WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji) 631 WRITE(numout,*) ' zqti_in : ', zqti_in (ji) * r1_rdtice 632 WRITE(numout,*) ' zqti_fin : ', zqti_fin(ji) * r1_rdtice 633 WRITE(numout,*) ' dh_i_bott: ', dh_i_bott(ji) 634 WRITE(numout,*) ' dh_i_surf: ', dh_i_surf(ji) 635 WRITE(numout,*) ' dh_snowice:', dh_snowice(ji) 636 WRITE(numout,*) ' icsuswi : ', icsuswi(ji) 637 WRITE(numout,*) ' icboswi : ', icboswi(ji) 638 WRITE(numout,*) ' snicswi : ', snicswi(ji) 639 ENDIF 640 END DO 641 ENDIF 641 642 642 643 !----------------------
Note: See TracChangeset
for help on using the changeset viewer.
