Changeset 4872 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
- Timestamp:
- 2014-11-18T18:03:00+01:00 (10 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 (modified) (45 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r4870 r4872 75 75 !! 76 76 !! ** Inputs / Ouputs : (global commons) 77 !! surface temperature : t_su_ b78 !! ice/snow temperatures : t_i_ b, t_s_b79 !! ice salinities : s_i_ b77 !! surface temperature : t_su_1d 78 !! ice/snow temperatures : t_i_1d, t_s_1d 79 !! ice salinities : s_i_1d 80 80 !! number of layers in the ice/snow: nlay_i, nlay_s 81 81 !! profile of the ice/snow layers : z_i, z_s 82 !! total ice/snow thickness : ht_i_ b, ht_s_b82 !! total ice/snow thickness : ht_i_1d, ht_s_1d 83 83 !! 84 84 !! ** External : … … 114 114 REAL(wp) :: zerritmax ! current maximal error on temperature 115 115 REAL(wp), POINTER, DIMENSION(:) :: ztfs ! ice melting point 116 REAL(wp), POINTER, DIMENSION(:) :: ztsu old! old surface temperature (before the iterative procedure )117 REAL(wp), POINTER, DIMENSION(:) :: ztsu oldit! surface temperature at previous iteration116 REAL(wp), POINTER, DIMENSION(:) :: ztsub ! old surface temperature (before the iterative procedure ) 117 REAL(wp), POINTER, DIMENSION(:) :: ztsubit ! surface temperature at previous iteration 118 118 REAL(wp), POINTER, DIMENSION(:) :: zh_i ! ice layer thickness 119 119 REAL(wp), POINTER, DIMENSION(:) :: zh_s ! snow layer thickness … … 129 129 REAL(wp), POINTER, DIMENSION(:,:) :: zradab_i ! Radiation absorbed in the ice 130 130 REAL(wp), POINTER, DIMENSION(:,:) :: zkappa_i ! Kappa factor in the ice 131 REAL(wp), POINTER, DIMENSION(:,:) :: zti old! Old temperature in the ice131 REAL(wp), POINTER, DIMENSION(:,:) :: ztib ! Old temperature in the ice 132 132 REAL(wp), POINTER, DIMENSION(:,:) :: zeta_i ! Eta factor in the ice 133 133 REAL(wp), POINTER, DIMENSION(:,:) :: ztitemp ! Temporary temperature in the ice to check the convergence … … 137 137 REAL(wp), POINTER, DIMENSION(:,:) :: zradab_s ! Radiation absorbed in the snow 138 138 REAL(wp), POINTER, DIMENSION(:,:) :: zkappa_s ! Kappa factor in the snow 139 REAL(wp), POINTER, DIMENSION(:,:) :: zeta_s 140 REAL(wp), POINTER, DIMENSION(:,:) :: ztstemp 141 REAL(wp), POINTER, DIMENSION(:,:) :: zts old! Temporary temperature in the snow142 REAL(wp), POINTER, DIMENSION(:,:) :: z_s 143 REAL(wp), POINTER, DIMENSION(:,:) :: zindterm ! Independent term144 REAL(wp), POINTER, DIMENSION(:,:) :: zindtbis ! temporary independent term139 REAL(wp), POINTER, DIMENSION(:,:) :: zeta_s ! Eta factor in the snow 140 REAL(wp), POINTER, DIMENSION(:,:) :: ztstemp ! Temporary temperature in the snow to check the convergence 141 REAL(wp), POINTER, DIMENSION(:,:) :: ztsb ! Temporary temperature in the snow 142 REAL(wp), POINTER, DIMENSION(:,:) :: z_s ! Vertical cotes of the layers in the snow 143 REAL(wp), POINTER, DIMENSION(:,:) :: zindterm ! Independent term 144 REAL(wp), POINTER, DIMENSION(:,:) :: zindtbis ! temporary independent term 145 145 REAL(wp), POINTER, DIMENSION(:,:) :: zdiagbis 146 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrid ! tridiagonal system terms146 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrid ! tridiagonal system terms 147 147 ! diag errors on heat 148 148 REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini … … 151 151 ! 152 152 CALL wrk_alloc( jpij, numeqmin, numeqmax, isnow ) 153 CALL wrk_alloc( jpij, ztfs, ztsu old, ztsuoldit, zh_i, zh_s, zfsw )153 CALL wrk_alloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw ) 154 154 CALL wrk_alloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu ) 155 CALL wrk_alloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, zti old, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0)156 CALL wrk_alloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, zts old, zeta_s, ztstemp, z_s, kjstart=0)155 CALL wrk_alloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0) 156 CALL wrk_alloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0) 157 157 CALL wrk_alloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis ) 158 158 CALL wrk_alloc( jpij, jkmax+2, 3, ztrid ) … … 163 163 zdq(:) = 0._wp ; zq_ini(:) = 0._wp 164 164 DO ji = kideb, kiut 165 zq_ini(ji) = ( SUM( q_i_ b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i ) + &166 & SUM( q_s_ b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ) )165 zq_ini(ji) = ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) / REAL( nlay_i ) + & 166 & SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) / REAL( nlay_s ) ) 167 167 END DO 168 168 … … 173 173 DO ji = kideb , kiut 174 174 ! is there snow or not 175 isnow(ji)= NINT( 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_ b(ji) ) ) )175 isnow(ji)= NINT( 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) ) ) 176 176 ! surface temperature of fusion 177 177 ztfs(ji) = REAL( isnow(ji) ) * rtt + REAL( 1 - isnow(ji) ) * rtt 178 178 ! layer thickness 179 zh_i(ji) = ht_i_ b(ji) / REAL( nlay_i )180 zh_s(ji) = ht_s_ b(ji) / REAL( nlay_s )179 zh_i(ji) = ht_i_1d(ji) / REAL( nlay_i ) 180 zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s ) 181 181 END DO 182 182 … … 190 190 DO jk = 1, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer 191 191 DO ji = kideb , kiut 192 z_s(ji,jk) = z_s(ji,jk-1) + ht_s_ b(ji) / REAL( nlay_s )192 z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) / REAL( nlay_s ) 193 193 END DO 194 194 END DO … … 196 196 DO jk = 1, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer 197 197 DO ji = kideb , kiut 198 z_i(ji,jk) = z_i(ji,jk-1) + ht_i_ b(ji) / REAL( nlay_i )198 z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) / REAL( nlay_i ) 199 199 END DO 200 200 END DO … … 217 217 DO ji = kideb , kiut 218 218 ! switches 219 isnow(ji) = NINT( 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_ b(ji) ) ) )219 isnow(ji) = NINT( 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) ) 220 220 ! hs > 0, isnow = 1 221 221 zhsu (ji) = hnzst ! threshold for the computation of i0 222 zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_ b(ji) / zhsu(ji) ) )222 zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) / zhsu(ji) ) ) 223 223 224 224 i0(ji) = REAL( 1 - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) ) … … 227 227 ! a function of the cloud cover 228 228 ! 229 !i0(ji) = (1.0-FLOAT(isnow(ji)))*3.0/(100*ht_s_ b(ji)+10.0)229 !i0(ji) = (1.0-FLOAT(isnow(ji)))*3.0/(100*ht_s_1d(ji)+10.0) 230 230 !formula used in Cice 231 231 END DO … … 272 272 273 273 DO ji = kideb, kiut ! Radiation transmitted below the ice 274 !!!ftr_ice_1d(ji) = ftr_ice_1d(ji) + iatte_1d(ji) * zradtr_i(ji,nlay_i) * a_i_ b(ji) / at_i_b(ji) ! clem modif274 !!!ftr_ice_1d(ji) = ftr_ice_1d(ji) + iatte_1d(ji) * zradtr_i(ji,nlay_i) * a_i_1d(ji) / at_i_1d(ji) ! clem modif 275 275 ftr_ice_1d(ji) = zradtr_i(ji,nlay_i) 276 276 END DO … … 282 282 ! 283 283 DO ji = kideb, kiut ! Old surface temperature 284 ztsu old (ji) = t_su_b(ji) ! temperature at the beg of iter pr.285 ztsu oldit(ji) = t_su_b(ji) ! temperature at the previous iter286 t_su_ b (ji) = MIN( t_su_b(ji), ztfs(ji) - ztsu_err ) ! necessary284 ztsub (ji) = t_su_1d(ji) ! temperature at the beg of iter pr. 285 ztsubit(ji) = t_su_1d(ji) ! temperature at the previous iter 286 t_su_1d (ji) = MIN( t_su_1d(ji), ztfs(ji) - ztsu_err ) ! necessary 287 287 zerrit (ji) = 1000._wp ! initial value of error 288 288 END DO … … 290 290 DO jk = 1, nlay_s ! Old snow temperature 291 291 DO ji = kideb , kiut 292 zts old(ji,jk) = t_s_b(ji,jk)292 ztsb(ji,jk) = t_s_1d(ji,jk) 293 293 END DO 294 294 END DO … … 296 296 DO jk = 1, nlay_i ! Old ice temperature 297 297 DO ji = kideb , kiut 298 zti old(ji,jk) = t_i_b(ji,jk)298 ztib(ji,jk) = t_i_1d(ji,jk) 299 299 END DO 300 300 END DO … … 313 313 IF( thcon_i_swi == 0 ) THEN ! Untersteiner (1964) formula 314 314 DO ji = kideb , kiut 315 ztcond_i(ji,0) = rcdic + zbeta*s_i_ b(ji,1) / MIN(-epsi10,t_i_b(ji,1)-rtt)315 ztcond_i(ji,0) = rcdic + zbeta*s_i_1d(ji,1) / MIN(-epsi10,t_i_1d(ji,1)-rtt) 316 316 ztcond_i(ji,0) = MAX(ztcond_i(ji,0),zkimin) 317 317 END DO 318 318 DO jk = 1, nlay_i-1 319 319 DO ji = kideb , kiut 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)320 ztcond_i(ji,jk) = rcdic + zbeta*( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) / & 321 MIN(-2.0_wp * epsi10, t_i_1d(ji,jk)+t_i_1d(ji,jk+1) - 2.0_wp * rtt) 322 322 ztcond_i(ji,jk) = MAX(ztcond_i(ji,jk),zkimin) 323 323 END DO … … 327 327 IF( thcon_i_swi == 1 ) THEN ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T 328 328 DO ji = kideb , kiut 329 ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_ b(ji,1) / MIN( -epsi10, t_i_b(ji,1)-rtt ) &330 & - 0.011_wp * ( t_i_ b(ji,1) - rtt )329 ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1)-rtt ) & 330 & - 0.011_wp * ( t_i_1d(ji,1) - rtt ) 331 331 ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) 332 332 END DO … … 334 334 DO ji = kideb , kiut 335 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 )336 & 0.090_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) & 337 & / MIN(-2.0_wp * epsi10, t_i_1d(ji,jk)+t_i_1d(ji,jk+1) - 2.0_wp * rtt) & 338 & - 0.0055_wp* ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0*rtt ) 339 339 ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin ) 340 340 END DO 341 341 END DO 342 342 DO ji = kideb , kiut 343 ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_ b(ji,nlay_i) / MIN(-epsi10,t_bo_b(ji)-rtt) &344 & - 0.011_wp * ( t_bo_ b(ji) - rtt )343 ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN(-epsi10,t_bo_1d(ji)-rtt) & 344 & - 0.011_wp * ( t_bo_1d(ji) - rtt ) 345 345 ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin ) 346 346 END DO … … 389 389 DO jk = 1, nlay_i 390 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)391 ztitemp(ji,jk) = t_i_1d(ji,jk) 392 zspeche_i(ji,jk) = cpic + zgamma*s_i_1d(ji,jk)/ & 393 MAX((t_i_1d(ji,jk)-rtt)*(ztib(ji,jk)-rtt),epsi10) 394 394 zeta_i(ji,jk) = rdt_ice / MAX(rhoic*zspeche_i(ji,jk)*zh_i(ji), & 395 395 epsi10) … … 399 399 DO jk = 1, nlay_s 400 400 DO ji = kideb , kiut 401 ztstemp(ji,jk) = t_s_ b(ji,jk)401 ztstemp(ji,jk) = t_s_1d(ji,jk) 402 402 zeta_s(ji,jk) = rdt_ice / MAX(rhosn*cpic*zh_s(ji),epsi10) 403 403 END DO … … 410 410 DO ji = kideb , kiut 411 411 ! update of the non solar flux according to the update in T_su 412 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_ b(ji) - ztsuoldit(ji) )412 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsubit(ji) ) 413 413 414 414 ! update incoming flux … … 448 448 zkappa_i(ji,jk)) 449 449 ztrid(ji,numeq,3) = - zeta_i(ji,jk)*zkappa_i(ji,jk) 450 zindterm(ji,numeq) = zti old(ji,jk) + zeta_i(ji,jk)* &450 zindterm(ji,numeq) = ztib(ji,jk) + zeta_i(ji,jk)* & 451 451 zradab_i(ji,jk) 452 452 END DO … … 460 460 + zkappa_i(ji,nlay_i-1) ) 461 461 ztrid(ji,numeq,3) = 0.0 462 zindterm(ji,numeq) = zti old(ji,nlay_i) + zeta_i(ji,nlay_i)* &462 zindterm(ji,numeq) = ztib(ji,nlay_i) + zeta_i(ji,nlay_i)* & 463 463 ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i)*zg1 & 464 * t_bo_ b(ji) )464 * t_bo_1d(ji) ) 465 465 ENDDO 466 466 467 467 468 468 DO ji = kideb , kiut 469 IF ( ht_s_ b(ji).gt.0.0 ) THEN469 IF ( ht_s_1d(ji).gt.0.0 ) THEN 470 470 ! 471 471 !------------------------------------------------------------------------------| … … 480 480 zkappa_s(ji,jk) ) 481 481 ztrid(ji,numeq,3) = - zeta_s(ji,jk)*zkappa_s(ji,jk) 482 zindterm(ji,numeq) = zts old(ji,jk) + zeta_s(ji,jk)* &482 zindterm(ji,numeq) = ztsb(ji,jk) + zeta_s(ji,jk)* & 483 483 zradab_s(ji,jk) 484 484 END DO … … 488 488 ztrid(ji,nlay_s+2,3) = 0.0 489 489 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1)* & 490 t_bo_ b(ji)490 t_bo_1d(ji) 491 491 ENDIF 492 492 493 IF ( t_su_ b(ji) .LT. rtt ) THEN493 IF ( t_su_1d(ji) .LT. rtt ) THEN 494 494 495 495 !------------------------------------------------------------------------------| … … 504 504 ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0) 505 505 ztrid(ji,1,3) = zg1s*zkappa_s(ji,0) 506 zindterm(ji,1) = dzf(ji)*t_su_ b(ji) - zf(ji)506 zindterm(ji,1) = dzf(ji)*t_su_1d(ji) - zf(ji) 507 507 508 508 !!first layer of snow equation … … 510 510 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1)*(zkappa_s(ji,1) + zkappa_s(ji,0)*zg1s) 511 511 ztrid(ji,2,3) = - zeta_s(ji,1)* zkappa_s(ji,1) 512 zindterm(ji,2) = zts old(ji,1) + zeta_s(ji,1)*zradab_s(ji,1)512 zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1)*zradab_s(ji,1) 513 513 514 514 ELSE … … 527 527 zkappa_s(ji,0) * zg1s ) 528 528 ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1) 529 zindterm(ji,2) = zts old(ji,1) + zeta_s(ji,1) * &529 zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) * & 530 530 ( zradab_s(ji,1) + & 531 zkappa_s(ji,0) * zg1s * t_su_ b(ji) )531 zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 532 532 ENDIF 533 533 ELSE … … 537 537 !------------------------------------------------------------------------------| 538 538 ! 539 IF (t_su_ b(ji) .LT. rtt) THEN539 IF (t_su_1d(ji) .LT. rtt) THEN 540 540 ! 541 541 !------------------------------------------------------------------------------| … … 551 551 ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*zg1 552 552 ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*zg1 553 zindterm(ji,numeqmin(ji)) = dzf(ji)*t_su_ b(ji) - zf(ji)553 zindterm(ji,numeqmin(ji)) = dzf(ji)*t_su_1d(ji) - zf(ji) 554 554 555 555 !!first layer of ice equation … … 558 558 + zkappa_i(ji,0) * zg1 ) 559 559 ztrid(ji,numeqmin(ji)+1,3) = - zeta_i(ji,1)*zkappa_i(ji,1) 560 zindterm(ji,numeqmin(ji)+1)= zti old(ji,1) + zeta_i(ji,1)*zradab_i(ji,1)560 zindterm(ji,numeqmin(ji)+1)= ztib(ji,1) + zeta_i(ji,1)*zradab_i(ji,1) 561 561 562 562 !!case of only one layer in the ice (surface & ice equations are altered) … … 571 571 ztrid(ji,numeqmin(ji)+1,3) = 0.0 572 572 573 zindterm(ji,numeqmin(ji)+1) = zti old(ji,1) + zeta_i(ji,1)* &574 ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_ b(ji) )573 zindterm(ji,numeqmin(ji)+1) = ztib(ji,1) + zeta_i(ji,1)* & 574 ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji) ) 575 575 ENDIF 576 576 … … 591 591 zg1) 592 592 ztrid(ji,numeqmin(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1) 593 zindterm(ji,numeqmin(ji)) = zti old(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + &594 zkappa_i(ji,0) * zg1 * t_su_ b(ji) )593 zindterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + & 594 zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 595 595 596 596 !!case of only one layer in the ice (surface & ice equations are altered) … … 600 600 zkappa_i(ji,1)) 601 601 ztrid(ji,numeqmin(ji),3) = 0.0 602 zindterm(ji,numeqmin(ji)) = zti old(ji,1) + zeta_i(ji,1)* &603 (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_ b(ji)) &604 + t_su_ b(ji)*zeta_i(ji,1)*zkappa_i(ji,0)*2.0602 zindterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1)* & 603 (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji)) & 604 + t_su_1d(ji)*zeta_i(ji,1)*zkappa_i(ji,0)*2.0 605 605 ENDIF 606 606 … … 642 642 DO ji = kideb , kiut 643 643 ! ice temperatures 644 t_i_ b(ji,nlay_i) = zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji))644 t_i_1d(ji,nlay_i) = zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji)) 645 645 END DO 646 646 … … 648 648 DO ji = kideb , kiut 649 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)650 t_i_1d(ji,jk) = (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* & 651 t_i_1d(ji,jk+1))/zdiagbis(ji,numeq) 652 652 END DO 653 653 END DO … … 655 655 DO ji = kideb , kiut 656 656 ! snow temperatures 657 IF (ht_s_ b(ji).GT.0._wp) &658 t_s_ b(ji,nlay_s) = (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) &659 * t_i_ b(ji,1))/zdiagbis(ji,nlay_s+1) &660 * MAX(0.0,SIGN(1.0,ht_s_ b(ji)))657 IF (ht_s_1d(ji).GT.0._wp) & 658 t_s_1d(ji,nlay_s) = (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) & 659 * t_i_1d(ji,1))/zdiagbis(ji,nlay_s+1) & 660 * MAX(0.0,SIGN(1.0,ht_s_1d(ji))) 661 661 662 662 ! surface temperature 663 isnow(ji) = NINT( 1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_ b(ji) ) ) )664 ztsu oldit(ji) = t_su_b(ji)665 IF( t_su_ b(ji) < ztfs(ji) ) &666 t_su_ b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_b(ji,1) &667 & + REAL( 1 - isnow(ji) )*t_i_ b(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))663 isnow(ji) = NINT( 1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_1d(ji) ) ) ) 664 ztsubit(ji) = t_su_1d(ji) 665 IF( t_su_1d(ji) < ztfs(ji) ) & 666 t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_1d(ji,1) & 667 & + REAL( 1 - isnow(ji) )*t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji)) 668 668 END DO 669 669 ! … … 675 675 ! zerrit(ji) is a measure of error, it has to be under maxer_i_thd 676 676 DO ji = kideb , kiut 677 t_su_ b(ji) = MAX( MIN( t_su_b(ji) , ztfs(ji) ) , 190._wp )678 zerrit(ji) = ABS( t_su_ b(ji) - ztsuoldit(ji) )677 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , ztfs(ji) ) , 190._wp ) 678 zerrit(ji) = ABS( t_su_1d(ji) - ztsubit(ji) ) 679 679 END DO 680 680 681 681 DO jk = 1, nlay_s 682 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)))683 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rtt ), 190._wp ) 684 zerrit(ji) = MAX(zerrit(ji),ABS(t_s_1d(ji,jk) - ztstemp(ji,jk))) 685 685 END DO 686 686 END DO … … 688 688 DO jk = 1, nlay_i 689 689 DO ji = kideb , kiut 690 ztmelt_i = -tmut * s_i_ b(ji,jk) + rtt691 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)))690 ztmelt_i = -tmut * s_i_1d(ji,jk) + rtt 691 t_i_1d(ji,jk) = MAX(MIN(t_i_1d(ji,jk),ztmelt_i), 190._wp) 692 zerrit(ji) = MAX(zerrit(ji),ABS(t_i_1d(ji,jk) - ztitemp(ji,jk))) 693 693 END DO 694 694 END DO … … 715 715 DO ji = kideb, kiut 716 716 ! forced mode only : update of latent heat fluxes (sublimation) (always >=0, upward flux) 717 IF( .NOT. lk_cpl) qla_ice_1d (ji) = MAX( 0._wp, qla_ice_1d (ji) + dqla_ice_1d(ji) * ( t_su_ b(ji) - ztsuold(ji) ) )717 IF( .NOT. lk_cpl) qla_ice_1d (ji) = MAX( 0._wp, qla_ice_1d (ji) + dqla_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) ) 718 718 ! ! surface ice conduction flux 719 isnow(ji) = NINT( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_ b(ji) ) ) )720 fc_su(ji) = - REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * (t_s_ b(ji,1) - t_su_b(ji)) &721 & - REAL( 1 - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_ b(ji,1) - t_su_b(ji))719 isnow(ji) = NINT( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) ) ) 720 fc_su(ji) = - REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji)) & 721 & - REAL( 1 - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_1d(ji,1) - t_su_1d(ji)) 722 722 ! ! bottom ice conduction flux 723 fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_ b(ji) - t_i_b(ji,nlay_i)) )723 fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) ) 724 724 END DO 725 725 … … 728 728 !----------------------------------------- 729 729 DO ji = kideb, kiut 730 IF( t_su_ b(ji) < rtt ) THEN ! case T_su < 0degC730 IF( t_su_1d(ji) < rtt ) THEN ! case T_su < 0degC 731 731 hfx_dif_1d(ji) = hfx_dif_1d(ji) + & 732 & ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_ b(ji)732 & ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 733 733 ELSE ! case T_su = 0degC 734 734 hfx_dif_1d(ji) = hfx_dif_1d(ji) + & 735 & ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_ b(ji)735 & ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 736 736 ENDIF 737 737 END DO … … 742 742 ! --- diag error on heat diffusion - PART 2 --- ! 743 743 DO ji = kideb, kiut 744 zdq(ji) = - zq_ini(ji) + ( SUM( q_i_ b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i ) + &745 & SUM( q_s_ b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ) )744 zdq(ji) = - zq_ini(ji) + ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) / REAL( nlay_i ) + & 745 & SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) / REAL( nlay_s ) ) 746 746 zhfx_err = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice ) 747 hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err * a_i_ b(ji)747 hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err * a_i_1d(ji) 748 748 ! --- correction of qns_ice and surface conduction flux --- ! 749 749 qns_ice_1d(ji) = qns_ice_1d(ji) - zhfx_err … … 751 751 ! --- Heat flux at the ice surface in W.m-2 --- ! 752 752 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 753 hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_ b(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) )753 hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_1d(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) ) 754 754 END DO 755 755 756 756 ! 757 757 CALL wrk_dealloc( jpij, numeqmin, numeqmax, isnow ) 758 CALL wrk_dealloc( jpij, ztfs, ztsu old, ztsuoldit, zh_i, zh_s, zfsw )758 CALL wrk_dealloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw ) 759 759 CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu ) 760 760 CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, & 761 & zti old, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 )762 CALL wrk_dealloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, zts old, zeta_s, ztstemp, z_s, kjstart = 0 )761 & ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 762 CALL wrk_dealloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 ) 763 763 CALL wrk_dealloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis ) 764 764 CALL wrk_dealloc( jpij, jkmax+2, 3, ztrid ) … … 783 783 DO jk = 1, nlay_i ! Sea ice energy of melting 784 784 DO ji = kideb, kiut 785 ztmelts = - tmut * s_i_ b(ji,jk) + rtt786 zindb = MAX( 0._wp , SIGN( 1._wp , -(t_i_ b(ji,jk) - rtt) - epsi10 ) )787 q_i_ b(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_b(ji,jk) ) &788 & + lfus * ( 1.0 - zindb * ( ztmelts-rtt ) / MIN( t_i_ b(ji,jk)-rtt, -epsi10 ) ) &785 ztmelts = - tmut * s_i_1d(ji,jk) + rtt 786 zindb = MAX( 0._wp , SIGN( 1._wp , -(t_i_1d(ji,jk) - rtt) - epsi10 ) ) 787 q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) ) & 788 & + lfus * ( 1.0 - zindb * ( ztmelts-rtt ) / MIN( t_i_1d(ji,jk)-rtt, -epsi10 ) ) & 789 789 & - rcp * ( ztmelts-rtt ) ) 790 790 END DO … … 792 792 DO jk = 1, nlay_s ! Snow energy of melting 793 793 DO ji = kideb, kiut 794 q_s_ b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )794 q_s_1d(ji,jk) = rhosn * ( cpic * ( rtt - t_s_1d(ji,jk) ) + lfus ) 795 795 END DO 796 796 END DO
Note: See TracChangeset
for help on using the changeset viewer.
