- Timestamp:
- 2014-11-27T17:13:38+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r4901 r4902 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 : … … 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 … … 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, zhfx_err … … 150 150 ! 151 151 CALL wrk_alloc( jpij, numeqmin, numeqmax, isnow ) 152 CALL wrk_alloc( jpij, ztfs, ztsu old, ztsuoldit, zh_i, zh_s, zfsw )152 CALL wrk_alloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw ) 153 153 CALL wrk_alloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu ) 154 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)155 CALL wrk_alloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, zts old, zeta_s, ztstemp, z_s, kjstart=0)156 CALL wrk_alloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis )157 CALL wrk_alloc( jpij, jkmax+2, 3, ztrid )154 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) 155 CALL wrk_alloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0) 156 CALL wrk_alloc( jpij, nlay_i+3, zindterm, zindtbis, zdiagbis ) 157 CALL wrk_alloc( jpij, nlay_i+3, 3, ztrid ) 158 158 159 159 CALL wrk_alloc( jpij, zdq, zq_ini, zhfx_err ) … … 162 162 zdq(:) = 0._wp ; zq_ini(:) = 0._wp 163 163 DO ji = kideb, kiut 164 zq_ini(ji) = ( SUM( q_i_ b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i ) + &165 & SUM( q_s_ b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ) )164 zq_ini(ji) = ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) / REAL( nlay_i ) + & 165 & SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) / REAL( nlay_s ) ) 166 166 END DO 167 167 … … 172 172 DO ji = kideb , kiut 173 173 ! is there snow or not 174 isnow(ji)= NINT( 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_ b(ji) ) ) )174 isnow(ji)= NINT( 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) ) ) 175 175 ! surface temperature of fusion 176 176 ztfs(ji) = REAL( isnow(ji) ) * rtt + REAL( 1 - isnow(ji) ) * rtt 177 177 ! layer thickness 178 zh_i(ji) = ht_i_ b(ji) / REAL( nlay_i )179 zh_s(ji) = ht_s_ b(ji) / REAL( nlay_s )178 zh_i(ji) = ht_i_1d(ji) / REAL( nlay_i ) 179 zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s ) 180 180 END DO 181 181 … … 187 187 z_i(:,0) = 0._wp ! vert. coord. of the up. lim. of the 1st ice layer 188 188 189 DO layer= 1, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer190 DO ji = kideb , kiut 191 z_s(ji, layer) = z_s(ji,layer-1) + ht_s_b(ji) / REAL( nlay_s )192 END DO 193 END DO 194 195 DO layer= 1, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer196 DO ji = kideb , kiut 197 z_i(ji, layer) = z_i(ji,layer-1) + ht_i_b(ji) / REAL( nlay_i )189 DO jk = 1, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer 190 DO ji = kideb , kiut 191 z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) / REAL( nlay_s ) 192 END DO 193 END DO 194 195 DO jk = 1, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer 196 DO ji = kideb , kiut 197 z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) / REAL( nlay_i ) 198 198 END DO 199 199 END DO … … 216 216 DO ji = kideb , kiut 217 217 ! switches 218 isnow(ji) = NINT( 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_ b(ji) ) ) )218 isnow(ji) = NINT( 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) ) 219 219 ! hs > 0, isnow = 1 220 220 zhsu (ji) = hnzst ! threshold for the computation of i0 221 zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_ b(ji) / zhsu(ji) ) )221 zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) / zhsu(ji) ) ) 222 222 223 223 i0(ji) = REAL( 1 - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) ) … … 226 226 ! a function of the cloud cover 227 227 ! 228 !i0(ji) = (1.0-FLOAT(isnow(ji)))*3.0/(100*ht_s_ b(ji)+10.0)228 !i0(ji) = (1.0-FLOAT(isnow(ji)))*3.0/(100*ht_s_1d(ji)+10.0) 229 229 !formula used in Cice 230 230 END DO … … 248 248 END DO 249 249 250 DO layer= 1, nlay_s ! Radiation through snow250 DO jk = 1, nlay_s ! Radiation through snow 251 251 DO ji = kideb, kiut 252 252 ! ! radiation transmitted below the layer-th snow layer 253 zradtr_s(ji, layer) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,layer) ) ) )253 zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,jk) ) ) ) 254 254 ! ! radiation absorbed by the layer-th snow layer 255 zradab_s(ji, layer) = zradtr_s(ji,layer-1) - zradtr_s(ji,layer)255 zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk) 256 256 END DO 257 257 END DO … … 261 261 END DO 262 262 263 DO layer= 1, nlay_i ! Radiation through ice263 DO jk = 1, nlay_i ! Radiation through ice 264 264 DO ji = kideb, kiut 265 265 ! ! radiation transmitted below the layer-th ice layer 266 zradtr_i(ji, layer) = zradtr_i(ji,0) * EXP( - kappa_i * ( MAX ( 0._wp , z_i(ji,layer) ) ) )266 zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - kappa_i * ( MAX ( 0._wp , z_i(ji,jk) ) ) ) 267 267 ! ! radiation absorbed by the layer-th ice layer 268 zradab_i(ji, layer) = zradtr_i(ji,layer-1) - zradtr_i(ji,layer)268 zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk) 269 269 END DO 270 270 END DO … … 280 280 ! 281 281 DO ji = kideb, kiut ! Old surface temperature 282 ztsu old (ji) = t_su_b(ji) ! temperature at the beg of iter pr.283 ztsu oldit(ji) = t_su_b(ji) ! temperature at the previous iter284 t_su_ b (ji) = MIN( t_su_b(ji), ztfs(ji) - ztsu_err ) ! necessary282 ztsub (ji) = t_su_1d(ji) ! temperature at the beg of iter pr. 283 ztsubit(ji) = t_su_1d(ji) ! temperature at the previous iter 284 t_su_1d (ji) = MIN( t_su_1d(ji), ztfs(ji) - ztsu_err ) ! necessary 285 285 zerrit (ji) = 1000._wp ! initial value of error 286 286 END DO 287 287 288 DO layer= 1, nlay_s ! Old snow temperature289 DO ji = kideb , kiut 290 zts old(ji,layer) = t_s_b(ji,layer)291 END DO 292 END DO 293 294 DO layer= 1, nlay_i ! Old ice temperature295 DO ji = kideb , kiut 296 zti old(ji,layer) = t_i_b(ji,layer)288 DO jk = 1, nlay_s ! Old snow temperature 289 DO ji = kideb , kiut 290 ztsb(ji,jk) = t_s_1d(ji,jk) 291 END DO 292 END DO 293 294 DO jk = 1, nlay_i ! Old ice temperature 295 DO ji = kideb , kiut 296 ztib(ji,jk) = t_i_1d(ji,jk) 297 297 END DO 298 298 END DO … … 311 311 IF( thcon_i_swi == 0 ) THEN ! Untersteiner (1964) formula 312 312 DO ji = kideb , kiut 313 ztcond_i(ji,0) = rcdic + zbeta*s_i_ b(ji,1) / MIN(-epsi10,t_i_b(ji,1)-rtt)313 ztcond_i(ji,0) = rcdic + zbeta*s_i_1d(ji,1) / MIN(-epsi10,t_i_1d(ji,1)-rtt) 314 314 ztcond_i(ji,0) = MAX(ztcond_i(ji,0),zkimin) 315 315 END DO 316 DO layer= 1, nlay_i-1316 DO jk = 1, nlay_i-1 317 317 DO ji = kideb , kiut 318 ztcond_i(ji, layer) = rcdic + zbeta*( s_i_b(ji,layer) + s_i_b(ji,layer+1) ) / &319 MIN(-2.0_wp * epsi10, t_i_ b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)320 ztcond_i(ji, layer) = MAX(ztcond_i(ji,layer),zkimin)318 ztcond_i(ji,jk) = rcdic + zbeta*( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) / & 319 MIN(-2.0_wp * epsi10, t_i_1d(ji,jk)+t_i_1d(ji,jk+1) - 2.0_wp * rtt) 320 ztcond_i(ji,jk) = MAX(ztcond_i(ji,jk),zkimin) 321 321 END DO 322 322 END DO … … 325 325 IF( thcon_i_swi == 1 ) THEN ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T 326 326 DO ji = kideb , kiut 327 ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_ b(ji,1) / MIN( -epsi10, t_i_b(ji,1)-rtt ) &328 & - 0.011_wp * ( t_i_ b(ji,1) - rtt )327 ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1)-rtt ) & 328 & - 0.011_wp * ( t_i_1d(ji,1) - rtt ) 329 329 ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) 330 330 END DO 331 DO layer= 1, nlay_i-1331 DO jk = 1, nlay_i-1 332 332 DO ji = kideb , kiut 333 ztcond_i(ji,layer) = rcdic + 0.090_wp * ( s_i_b(ji,layer) + s_i_b(ji,layer+1) ) & 334 & / MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt) & 335 & - 0.0055_wp* ( t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2.0*rtt ) 336 ztcond_i(ji,layer) = MAX( ztcond_i(ji,layer), zkimin ) 333 ztcond_i(ji,jk) = rcdic + & 334 & 0.090_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) & 335 & / MIN(-2.0_wp * epsi10, t_i_1d(ji,jk)+t_i_1d(ji,jk+1) - 2.0_wp * rtt) & 336 & - 0.0055_wp* ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0*rtt ) 337 ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin ) 337 338 END DO 338 339 END DO 339 340 DO ji = kideb , kiut 340 ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_ b(ji,nlay_i) / MIN(-epsi10,t_bo_b(ji)-rtt) &341 & - 0.011_wp * ( t_bo_ b(ji) - rtt )341 ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN(-epsi10,t_bo_1d(ji)-rtt) & 342 & - 0.011_wp * ( t_bo_1d(ji) - rtt ) 342 343 ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin ) 343 344 END DO … … 355 356 END DO 356 357 357 DO layer= 1, nlay_s-1358 DO ji = kideb , kiut 359 zkappa_s(ji, layer) = 2.0 * rcdsn / &358 DO jk = 1, nlay_s-1 359 DO ji = kideb , kiut 360 zkappa_s(ji,jk) = 2.0 * rcdsn / & 360 361 MAX(epsi10,2.0*zh_s(ji)) 361 362 END DO 362 363 END DO 363 364 364 DO layer= 1, nlay_i-1365 DO jk = 1, nlay_i-1 365 366 DO ji = kideb , kiut 366 367 !-- Ice kappa factors 367 zkappa_i(ji, layer) = 2.0*ztcond_i(ji,layer)/ &368 zkappa_i(ji,jk) = 2.0*ztcond_i(ji,jk)/ & 368 369 MAX(epsi10,2.0*zh_i(ji)) 369 370 END DO … … 384 385 !------------------------------------------------------------------------------| 385 386 ! 386 DO layer= 1, nlay_i387 DO ji = kideb , kiut 388 ztitemp(ji, layer) = t_i_b(ji,layer)389 zspeche_i(ji, layer) = cpic + zgamma*s_i_b(ji,layer)/ &390 MAX((t_i_ b(ji,layer)-rtt)*(ztiold(ji,layer)-rtt),epsi10)391 zeta_i(ji, layer) = rdt_ice / MAX(rhoic*zspeche_i(ji,layer)*zh_i(ji), &387 DO jk = 1, nlay_i 388 DO ji = kideb , kiut 389 ztitemp(ji,jk) = t_i_1d(ji,jk) 390 zspeche_i(ji,jk) = cpic + zgamma*s_i_1d(ji,jk)/ & 391 MAX((t_i_1d(ji,jk)-rtt)*(ztib(ji,jk)-rtt),epsi10) 392 zeta_i(ji,jk) = rdt_ice / MAX(rhoic*zspeche_i(ji,jk)*zh_i(ji), & 392 393 epsi10) 393 394 END DO 394 395 END DO 395 396 396 DO layer= 1, nlay_s397 DO ji = kideb , kiut 398 ztstemp(ji, layer) = t_s_b(ji,layer)399 zeta_s(ji, layer) = rdt_ice / MAX(rhosn*cpic*zh_s(ji),epsi10)397 DO jk = 1, nlay_s 398 DO ji = kideb , kiut 399 ztstemp(ji,jk) = t_s_1d(ji,jk) 400 zeta_s(ji,jk) = rdt_ice / MAX(rhosn*cpic*zh_s(ji),epsi10) 400 401 END DO 401 402 END DO … … 408 409 DO ji = kideb , kiut 409 410 ! update of the non solar flux according to the update in T_su 410 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_ b(ji) - ztsuoldit(ji) )411 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsubit(ji) ) 411 412 END DO 412 413 ENDIF … … 432 433 !!ice interior terms (top equation has the same form as the others) 433 434 434 DO numeq=1, jkmax+2435 DO numeq=1,nlay_i+3 435 436 DO ji = kideb , kiut 436 437 ztrid(ji,numeq,1) = 0. … … 445 446 DO numeq = nlay_s + 2, nlay_s + nlay_i 446 447 DO ji = kideb , kiut 447 layer= numeq - nlay_s - 1448 ztrid(ji,numeq,1) = - zeta_i(ji, layer)*zkappa_i(ji,layer-1)449 ztrid(ji,numeq,2) = 1.0 + zeta_i(ji, layer)*(zkappa_i(ji,layer-1) + &450 zkappa_i(ji, layer))451 ztrid(ji,numeq,3) = - zeta_i(ji, layer)*zkappa_i(ji,layer)452 zindterm(ji,numeq) = zti old(ji,layer) + zeta_i(ji,layer)* &453 zradab_i(ji, layer)448 jk = numeq - nlay_s - 1 449 ztrid(ji,numeq,1) = - zeta_i(ji,jk)*zkappa_i(ji,jk-1) 450 ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,jk)*(zkappa_i(ji,jk-1) + & 451 zkappa_i(ji,jk)) 452 ztrid(ji,numeq,3) = - zeta_i(ji,jk)*zkappa_i(ji,jk) 453 zindterm(ji,numeq) = ztib(ji,jk) + zeta_i(ji,jk)* & 454 zradab_i(ji,jk) 454 455 END DO 455 456 ENDDO … … 462 463 + zkappa_i(ji,nlay_i-1) ) 463 464 ztrid(ji,numeq,3) = 0.0 464 zindterm(ji,numeq) = zti old(ji,nlay_i) + zeta_i(ji,nlay_i)* &465 zindterm(ji,numeq) = ztib(ji,nlay_i) + zeta_i(ji,nlay_i)* & 465 466 ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i)*zg1 & 466 * t_bo_ b(ji) )467 * t_bo_1d(ji) ) 467 468 ENDDO 468 469 469 470 470 471 DO ji = kideb , kiut 471 IF ( ht_s_ b(ji).gt.0.0 ) THEN472 IF ( ht_s_1d(ji).gt.0.0 ) THEN 472 473 ! 473 474 !------------------------------------------------------------------------------| … … 477 478 !!snow interior terms (bottom equation has the same form as the others) 478 479 DO numeq = 3, nlay_s + 1 479 layer= numeq - 1480 ztrid(ji,numeq,1) = - zeta_s(ji, layer)*zkappa_s(ji,layer-1)481 ztrid(ji,numeq,2) = 1.0 + zeta_s(ji, layer)*( zkappa_s(ji,layer-1) + &482 zkappa_s(ji, layer) )483 ztrid(ji,numeq,3) = - zeta_s(ji, layer)*zkappa_s(ji,layer)484 zindterm(ji,numeq) = zts old(ji,layer) + zeta_s(ji,layer)* &485 zradab_s(ji, layer)480 jk = numeq - 1 481 ztrid(ji,numeq,1) = - zeta_s(ji,jk)*zkappa_s(ji,jk-1) 482 ztrid(ji,numeq,2) = 1.0 + zeta_s(ji,jk)*( zkappa_s(ji,jk-1) + & 483 zkappa_s(ji,jk) ) 484 ztrid(ji,numeq,3) = - zeta_s(ji,jk)*zkappa_s(ji,jk) 485 zindterm(ji,numeq) = ztsb(ji,jk) + zeta_s(ji,jk)* & 486 zradab_s(ji,jk) 486 487 END DO 487 488 … … 490 491 ztrid(ji,nlay_s+2,3) = 0.0 491 492 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1)* & 492 t_bo_ b(ji)493 t_bo_1d(ji) 493 494 ENDIF 494 495 495 IF ( t_su_ b(ji) .LT. rtt ) THEN496 IF ( t_su_1d(ji) .LT. rtt ) THEN 496 497 497 498 !------------------------------------------------------------------------------| … … 506 507 ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0) 507 508 ztrid(ji,1,3) = zg1s*zkappa_s(ji,0) 508 zindterm(ji,1) = dzf(ji)*t_su_ b(ji) - zf(ji)509 zindterm(ji,1) = dzf(ji)*t_su_1d(ji) - zf(ji) 509 510 510 511 !!first layer of snow equation … … 512 513 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1)*(zkappa_s(ji,1) + zkappa_s(ji,0)*zg1s) 513 514 ztrid(ji,2,3) = - zeta_s(ji,1)* zkappa_s(ji,1) 514 zindterm(ji,2) = zts old(ji,1) + zeta_s(ji,1)*zradab_s(ji,1)515 zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1)*zradab_s(ji,1) 515 516 516 517 ELSE … … 529 530 zkappa_s(ji,0) * zg1s ) 530 531 ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1) 531 zindterm(ji,2) = zts old(ji,1) + zeta_s(ji,1) * &532 zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) * & 532 533 ( zradab_s(ji,1) + & 533 zkappa_s(ji,0) * zg1s * t_su_ b(ji) )534 zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 534 535 ENDIF 535 536 ELSE … … 539 540 !------------------------------------------------------------------------------| 540 541 ! 541 IF (t_su_ b(ji) .LT. rtt) THEN542 IF (t_su_1d(ji) .LT. rtt) THEN 542 543 ! 543 544 !------------------------------------------------------------------------------| … … 553 554 ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*zg1 554 555 ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*zg1 555 zindterm(ji,numeqmin(ji)) = dzf(ji)*t_su_ b(ji) - zf(ji)556 zindterm(ji,numeqmin(ji)) = dzf(ji)*t_su_1d(ji) - zf(ji) 556 557 557 558 !!first layer of ice equation … … 560 561 + zkappa_i(ji,0) * zg1 ) 561 562 ztrid(ji,numeqmin(ji)+1,3) = - zeta_i(ji,1)*zkappa_i(ji,1) 562 zindterm(ji,numeqmin(ji)+1)= zti old(ji,1) + zeta_i(ji,1)*zradab_i(ji,1)563 zindterm(ji,numeqmin(ji)+1)= ztib(ji,1) + zeta_i(ji,1)*zradab_i(ji,1) 563 564 564 565 !!case of only one layer in the ice (surface & ice equations are altered) … … 573 574 ztrid(ji,numeqmin(ji)+1,3) = 0.0 574 575 575 zindterm(ji,numeqmin(ji)+1) = zti old(ji,1) + zeta_i(ji,1)* &576 ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_ b(ji) )576 zindterm(ji,numeqmin(ji)+1) = ztib(ji,1) + zeta_i(ji,1)* & 577 ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji) ) 577 578 ENDIF 578 579 … … 593 594 zg1) 594 595 ztrid(ji,numeqmin(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1) 595 zindterm(ji,numeqmin(ji)) = zti old(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + &596 zkappa_i(ji,0) * zg1 * t_su_ b(ji) )596 zindterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + & 597 zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 597 598 598 599 !!case of only one layer in the ice (surface & ice equations are altered) … … 602 603 zkappa_i(ji,1)) 603 604 ztrid(ji,numeqmin(ji),3) = 0.0 604 zindterm(ji,numeqmin(ji)) = zti old(ji,1) + zeta_i(ji,1)* &605 (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_ b(ji)) &606 + t_su_ b(ji)*zeta_i(ji,1)*zkappa_i(ji,0)*2.0605 zindterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1)* & 606 (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji)) & 607 + t_su_1d(ji)*zeta_i(ji,1)*zkappa_i(ji,0)*2.0 607 608 ENDIF 608 609 … … 623 624 624 625 maxnumeqmax = 0 625 minnumeqmin = jkmax+4626 minnumeqmin = nlay_i+5 626 627 627 628 DO ji = kideb , kiut … … 632 633 END DO 633 634 634 DO layer= minnumeqmin+1, maxnumeqmax635 DO ji = kideb , kiut 636 numeq = min(max(numeqmin(ji)+1, layer),numeqmax(ji))635 DO jk = minnumeqmin+1, maxnumeqmax 636 DO ji = kideb , kiut 637 numeq = min(max(numeqmin(ji)+1,jk),numeqmax(ji)) 637 638 zdiagbis(ji,numeq) = ztrid(ji,numeq,2) - ztrid(ji,numeq,1)* & 638 639 ztrid(ji,numeq-1,3)/zdiagbis(ji,numeq-1) … … 644 645 DO ji = kideb , kiut 645 646 ! ice temperatures 646 t_i_ b(ji,nlay_i) = zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji))647 t_i_1d(ji,nlay_i) = zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji)) 647 648 END DO 648 649 649 650 DO numeq = nlay_i + nlay_s + 1, nlay_s + 2, -1 650 651 DO ji = kideb , kiut 651 layer= numeq - nlay_s - 1652 t_i_ b(ji,layer) = (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* &653 t_i_ b(ji,layer+1))/zdiagbis(ji,numeq)652 jk = numeq - nlay_s - 1 653 t_i_1d(ji,jk) = (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* & 654 t_i_1d(ji,jk+1))/zdiagbis(ji,numeq) 654 655 END DO 655 656 END DO … … 657 658 DO ji = kideb , kiut 658 659 ! snow temperatures 659 IF (ht_s_ b(ji).GT.0._wp) &660 t_s_ b(ji,nlay_s) = (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) &661 * t_i_ b(ji,1))/zdiagbis(ji,nlay_s+1) &662 * MAX(0.0,SIGN(1.0,ht_s_ b(ji)))660 IF (ht_s_1d(ji).GT.0._wp) & 661 t_s_1d(ji,nlay_s) = (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) & 662 * t_i_1d(ji,1))/zdiagbis(ji,nlay_s+1) & 663 * MAX(0.0,SIGN(1.0,ht_s_1d(ji))) 663 664 664 665 ! surface temperature 665 isnow(ji) = NINT( 1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_ b(ji) ) ) )666 ztsu oldit(ji) = t_su_b(ji)667 IF( t_su_ b(ji) < ztfs(ji) ) &668 t_su_ b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_b(ji,1) &669 & + REAL( 1 - isnow(ji) )*t_i_ b(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))666 isnow(ji) = NINT( 1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_1d(ji) ) ) ) 667 ztsubit(ji) = t_su_1d(ji) 668 IF( t_su_1d(ji) < ztfs(ji) ) & 669 t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_1d(ji,1) & 670 & + REAL( 1 - isnow(ji) )*t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji)) 670 671 END DO 671 672 ! … … 677 678 ! zerrit(ji) is a measure of error, it has to be under maxer_i_thd 678 679 DO ji = kideb , kiut 679 t_su_ b(ji) = MAX( MIN( t_su_b(ji) , ztfs(ji) ) , 190._wp )680 zerrit(ji) = ABS( t_su_ b(ji) - ztsuoldit(ji) )681 END DO 682 683 DO layer= 1, nlay_s684 DO ji = kideb , kiut 685 t_s_ b(ji,layer) = MAX( MIN( t_s_b(ji,layer), rtt ), 190._wp )686 zerrit(ji) = MAX(zerrit(ji),ABS(t_s_ b(ji,layer) - ztstemp(ji,layer)))687 END DO 688 END DO 689 690 DO layer= 1, nlay_i691 DO ji = kideb , kiut 692 ztmelt_i = -tmut * s_i_ b(ji,layer) + rtt693 t_i_ b(ji,layer) = MAX(MIN(t_i_b(ji,layer),ztmelt_i), 190._wp)694 zerrit(ji) = MAX(zerrit(ji),ABS(t_i_ b(ji,layer) - ztitemp(ji,layer)))680 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , ztfs(ji) ) , 190._wp ) 681 zerrit(ji) = ABS( t_su_1d(ji) - ztsubit(ji) ) 682 END DO 683 684 DO jk = 1, nlay_s 685 DO ji = kideb , kiut 686 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rtt ), 190._wp ) 687 zerrit(ji) = MAX(zerrit(ji),ABS(t_s_1d(ji,jk) - ztstemp(ji,jk))) 688 END DO 689 END DO 690 691 DO jk = 1, nlay_i 692 DO ji = kideb , kiut 693 ztmelt_i = -tmut * s_i_1d(ji,jk) + rtt 694 t_i_1d(ji,jk) = MAX(MIN(t_i_1d(ji,jk),ztmelt_i), 190._wp) 695 zerrit(ji) = MAX(zerrit(ji),ABS(t_i_1d(ji,jk) - ztitemp(ji,jk))) 695 696 END DO 696 697 END DO … … 717 718 DO ji = kideb, kiut 718 719 ! forced mode only : update of latent heat fluxes (sublimation) (always >=0, upward flux) 719 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) ) )720 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) ) ) 720 721 ! ! surface ice conduction flux 721 isnow(ji) = NINT( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_ b(ji) ) ) )722 fc_su(ji) = - REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * (t_s_ b(ji,1) - t_su_b(ji)) &723 & - REAL( 1 - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_ b(ji,1) - t_su_b(ji))722 isnow(ji) = NINT( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) ) ) 723 fc_su(ji) = - REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji)) & 724 & - REAL( 1 - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_1d(ji,1) - t_su_1d(ji)) 724 725 ! ! bottom ice conduction flux 725 fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_ b(ji) - t_i_b(ji,nlay_i)) )726 fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) ) 726 727 END DO 727 728 … … 730 731 !----------------------------------------- 731 732 DO ji = kideb, kiut 732 IF( t_su_b(ji) < rtt ) THEN ! case T_su < 0degC 733 hfx_dif_1d(ji) = hfx_dif_1d(ji) + ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji) 733 IF( t_su_1d(ji) < rtt ) THEN ! case T_su < 0degC 734 hfx_dif_1d(ji) = hfx_dif_1d(ji) + & 735 & ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 734 736 ELSE ! case T_su = 0degC 735 hfx_dif_1d(ji) = hfx_dif_1d(ji) + ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji) 737 hfx_dif_1d(ji) = hfx_dif_1d(ji) + & 738 & ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 736 739 ENDIF 737 740 END DO … … 742 745 ! --- diag conservation imbalance on heat diffusion - PART 2 --- ! 743 746 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 ) )747 zdq(ji) = - zq_ini(ji) + ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) / REAL( nlay_i ) + & 748 & SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) / REAL( nlay_s ) ) 746 749 zhfx_err(ji) = ( 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(ji) * a_i_ b(ji)750 hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_1d(ji) 748 751 END DO 749 752 … … 767 770 DO ji = kideb, kiut 768 771 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 769 hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_ b(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) )772 hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_1d(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) ) 770 773 END DO 771 774 772 775 ! 773 776 CALL wrk_dealloc( jpij, numeqmin, numeqmax, isnow ) 774 CALL wrk_dealloc( jpij, ztfs, ztsu old, ztsuoldit, zh_i, zh_s, zfsw )777 CALL wrk_dealloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw ) 775 778 CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu ) 776 CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztiold, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 777 CALL wrk_dealloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsold, zeta_s, ztstemp, z_s, kjstart = 0 ) 778 CALL wrk_dealloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis ) 779 CALL wrk_dealloc( jpij, jkmax+2, 3, ztrid ) 779 CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, & 780 & ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 781 CALL wrk_dealloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 ) 782 CALL wrk_dealloc( jpij, nlay_i+3, zindterm, zindtbis, zdiagbis ) 783 CALL wrk_dealloc( jpij, nlay_i+3, 3, ztrid ) 780 784 CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx_err ) 781 785 … … 798 802 DO jk = 1, nlay_i ! Sea ice energy of melting 799 803 DO ji = kideb, kiut 800 ztmelts = - tmut * s_i_ b(ji,jk) + rtt801 zindb = MAX( 0._wp , SIGN( 1._wp , -(t_i_ b(ji,jk) - rtt) - epsi10 ) )802 q_i_ b(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_b(ji,jk) ) &803 & + lfus * ( 1.0 - zindb * ( ztmelts-rtt ) / MIN( t_i_ b(ji,jk)-rtt, -epsi10 ) ) &804 ztmelts = - tmut * s_i_1d(ji,jk) + rtt 805 zindb = MAX( 0._wp , SIGN( 1._wp , -(t_i_1d(ji,jk) - rtt) - epsi10 ) ) 806 q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) ) & 807 & + lfus * ( 1.0 - zindb * ( ztmelts-rtt ) / MIN( t_i_1d(ji,jk)-rtt, -epsi10 ) ) & 804 808 & - rcp * ( ztmelts-rtt ) ) 805 809 END DO … … 807 811 DO jk = 1, nlay_s ! Snow energy of melting 808 812 DO ji = kideb, kiut 809 q_s_ b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )813 q_s_1d(ji,jk) = rhosn * ( cpic * ( rtt - t_s_1d(ji,jk) ) + lfus ) 810 814 END DO 811 815 END DO
Note: See TracChangeset
for help on using the changeset viewer.
