- Timestamp:
- 2017-09-27T12:09:10+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90
r8564 r8565 144 144 145 145 ! --- diag error on heat diffusion - PART 1 --- ! 146 DO ji = 1, n idx146 DO ji = 1, npti 147 147 zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i + & 148 148 & SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) … … 152 152 ! 1) Initialization 153 153 !------------------ 154 DO ji = 1, n idx154 DO ji = 1, npti 155 155 isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - h_s_1d(ji) ) ) ! is there snow or not 156 156 ! layer thickness … … 159 159 END DO 160 160 ! 161 WHERE( zh_i(1:n idx) >= epsi10 ) ; z1_h_i(1:nidx) = 1._wp / zh_i(1:nidx)162 ELSEWHERE ; z1_h_i(1:n idx) = 0._wp161 WHERE( zh_i(1:npti) >= epsi10 ) ; z1_h_i(1:npti) = 1._wp / zh_i(1:npti) 162 ELSEWHERE ; z1_h_i(1:npti) = 0._wp 163 163 END WHERE 164 164 165 WHERE( zh_s(1:n idx) >= epsi10 ) ; z1_h_s(1:nidx) = 1._wp / zh_s(1:nidx)166 ELSEWHERE ; z1_h_s(1:n idx) = 0._wp165 WHERE( zh_s(1:npti) >= epsi10 ) ; z1_h_s(1:npti) = 1._wp / zh_s(1:npti) 166 ELSEWHERE ; z1_h_s(1:npti) = 0._wp 167 167 END WHERE 168 168 ! 169 169 ! temperatures 170 ztsub (1:n idx) = t_su_1d(1:nidx) ! temperature at the previous iteration171 ztsold (1:n idx,:) = t_s_1d(1:nidx,:) ! Old snow temperature172 ztiold (1:n idx,:) = t_i_1d(1:nidx,:) ! Old ice temperature173 t_su_1d(1:n idx) = MIN( t_su_1d(1:nidx), rt0 - ztsu_err ) ! necessary170 ztsub (1:npti) = t_su_1d(1:npti) ! temperature at the previous iteration 171 ztsold (1:npti,:) = t_s_1d(1:npti,:) ! Old snow temperature 172 ztiold (1:npti,:) = t_i_1d(1:npti,:) ! Old ice temperature 173 t_su_1d(1:npti) = MIN( t_su_1d(1:npti), rt0 - ztsu_err ) ! necessary 174 174 ! 175 175 !------------- … … 177 177 !------------- 178 178 z1_hsu = 1._wp / 0.1_wp ! threshold for the computation of i0 179 DO ji = 1, n idx179 DO ji = 1, npti 180 180 ! --- Computation of i0 --- ! 181 181 ! i0 describes the fraction of solar radiation which does not contribute … … 199 199 200 200 ! --- Transmission/absorption of solar radiation in the ice --- ! 201 zradtr_s(1:n idx,0) = zftrice(1:nidx)201 zradtr_s(1:npti,0) = zftrice(1:npti) 202 202 DO jk = 1, nlay_s 203 DO ji = 1, n idx203 DO ji = 1, npti 204 204 ! ! radiation transmitted below the layer-th snow layer 205 205 zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * zh_s(ji) * REAL(jk) ) … … 209 209 END DO 210 210 211 zradtr_i(1:n idx,0) = zradtr_s(1:nidx,nlay_s) * isnow(1:nidx) + zftrice(1:nidx) * ( 1._wp - isnow(1:nidx) )211 zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + zftrice(1:npti) * ( 1._wp - isnow(1:npti) ) 212 212 DO jk = 1, nlay_i 213 DO ji = 1, n idx213 DO ji = 1, npti 214 214 ! ! radiation transmitted below the layer-th ice layer 215 215 zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) ) … … 219 219 END DO 220 220 221 ftr_ice_1d(1:n idx) = zradtr_i(1:nidx,nlay_i) ! record radiation transmitted below the ice221 ftr_ice_1d(1:npti) = zradtr_i(1:npti,nlay_i) ! record radiation transmitted below the ice 222 222 ! 223 223 iconv = 0 ! number of iterations … … 228 228 iconv = iconv + 1 229 229 ! 230 ztib(1:n idx,:) = t_i_1d(1:nidx,:)231 ztsb(1:n idx,:) = t_s_1d(1:nidx,:)230 ztib(1:npti,:) = t_i_1d(1:npti,:) 231 ztsb(1:npti,:) = t_s_1d(1:npti,:) 232 232 ! 233 233 !-------------------------------- … … 236 236 IF( ln_cndi_U64 ) THEN !-- Untersteiner (1964) formula: k = k0 + beta.S/T 237 237 ! 238 DO ji = 1, n idx238 DO ji = 1, npti 239 239 ztcond_i(ji,0) = rcdic + zbeta * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) 240 240 ztcond_i(ji,nlay_i) = rcdic + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) 241 241 END DO 242 242 DO jk = 1, nlay_i-1 243 DO ji = 1, n idx243 DO ji = 1, npti 244 244 ztcond_i(ji,jk) = rcdic + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / & 245 245 & MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) … … 249 249 ELSEIF( ln_cndi_P07 ) THEN !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T 250 250 ! 251 DO ji = 1, n idx251 DO ji = 1, npti 252 252 ztcond_i(ji,0) = rcdic + 0.09_wp * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) & 253 253 & - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) … … 256 256 END DO 257 257 DO jk = 1, nlay_i-1 258 DO ji = 1, n idx258 DO ji = 1, npti 259 259 ztcond_i(ji,jk) = rcdic + 0.09_wp * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / & 260 260 & MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) & … … 264 264 ! 265 265 ENDIF 266 ztcond_i(1:n idx,:) = MAX( zkimin, ztcond_i(1:nidx,:) )266 ztcond_i(1:npti,:) = MAX( zkimin, ztcond_i(1:npti,:) ) 267 267 ! 268 268 !--- G(he) : enhancement of thermal conductivity in mono-category case … … 270 270 ! Used in mono-category case only to simulate an ITD implicitly 271 271 ! Fichefet and Morales Maqueda, JGR 1997 272 zghe(1:n idx) = 1._wp272 zghe(1:npti) = 1._wp 273 273 ! 274 274 SELECT CASE ( nn_monocat ) … … 277 277 278 278 zepsilon = 0.1_wp 279 DO ji = 1, n idx279 DO ji = 1, npti 280 280 zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp ) ! Mean sea ice thermal conductivity 281 281 zhe = ( rn_cnd_s * h_i_1d(ji) + zcnd_i * h_s_1d(ji) ) / ( rn_cnd_s + zcnd_i ) ! Effective thickness he (zhe) … … 292 292 !--- Snow 293 293 DO jk = 0, nlay_s-1 294 DO ji = 1, n idx294 DO ji = 1, npti 295 295 zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji) 296 296 END DO 297 297 END DO 298 DO ji = 1, n idx! Snow-ice interface298 DO ji = 1, npti ! Snow-ice interface 299 299 zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) 300 300 IF( zfac > epsi10 ) THEN … … 307 307 !--- Ice 308 308 DO jk = 0, nlay_i 309 DO ji = 1, n idx309 DO ji = 1, npti 310 310 zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji) 311 311 END DO 312 312 END DO 313 DO ji = 1, n idx! Snow-ice interface313 DO ji = 1, npti ! Snow-ice interface 314 314 zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 315 315 END DO … … 319 319 !-------------------------------------- 320 320 DO jk = 1, nlay_i 321 DO ji = 1, n idx321 DO ji = 1, npti 322 322 zcpi = cpic + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 ) 323 323 zeta_i(ji,jk) = rdt_ice * r1_rhoic * z1_h_i(ji) / MAX( epsi10, zcpi ) … … 326 326 327 327 DO jk = 1, nlay_s 328 DO ji = 1, n idx328 DO ji = 1, npti 329 329 zeta_s(ji,jk) = rdt_ice * r1_rhosn * r1_cpic * z1_h_s(ji) 330 330 END DO … … 335 335 !---------------------------- 336 336 IF ( ln_dqns_i ) THEN 337 DO ji = 1, n idx337 DO ji = 1, npti 338 338 ! update of the non solar flux according to the update in T_su 339 339 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) … … 341 341 ENDIF 342 342 343 DO ji = 1, n idx343 DO ji = 1, npti 344 344 zfnet(ji) = zfsw(ji) + qns_ice_1d(ji) ! incoming = net absorbed solar radiation + non solar total flux (LWup, LWdw, SH, LH) 345 345 END DO … … 354 354 355 355 !!ice interior terms (top equation has the same form as the others) 356 ztrid (1:n idx,:,:) = 0._wp357 zindterm(1:n idx,:) = 0._wp358 zindtbis(1:n idx,:) = 0._wp359 zdiagbis(1:n idx,:) = 0._wp356 ztrid (1:npti,:,:) = 0._wp 357 zindterm(1:npti,:) = 0._wp 358 zindtbis(1:npti,:) = 0._wp 359 zdiagbis(1:npti,:) = 0._wp 360 360 361 361 DO jm = nlay_s + 2, nlay_s + nlay_i 362 DO ji = 1, n idx362 DO ji = 1, npti 363 363 jk = jm - nlay_s - 1 364 364 ztrid(ji,jm,1) = - zeta_i(ji,jk) * zkappa_i(ji,jk-1) … … 370 370 371 371 jm = nlay_s + nlay_i + 1 372 DO ji = 1, n idx372 DO ji = 1, npti 373 373 !!ice bottom term 374 374 ztrid(ji,jm,1) = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1) … … 380 380 381 381 382 DO ji = 1, n idx382 DO ji = 1, npti 383 383 ! !---------------------! 384 384 IF ( h_s_1d(ji) > 0.0 ) THEN ! snow-covered cells ! … … 497 497 jm_maxt = 0 498 498 jm_mint = nlay_i+5 499 DO ji = 1, n idx499 DO ji = 1, npti 500 500 jm_mint = MIN(jm_min(ji),jm_mint) 501 501 jm_maxt = MAX(jm_max(ji),jm_maxt) … … 503 503 504 504 DO jk = jm_mint+1, jm_maxt 505 DO ji = 1, n idx505 DO ji = 1, npti 506 506 jm = min(max(jm_min(ji)+1,jk),jm_max(ji)) 507 507 zdiagbis(ji,jm) = ztrid(ji,jm,2) - ztrid(ji,jm,1) * ztrid(ji,jm-1,3) / zdiagbis(ji,jm-1) … … 510 510 END DO 511 511 512 DO ji = 1, n idx512 DO ji = 1, npti 513 513 ! ice temperatures 514 514 t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) … … 516 516 517 517 DO jm = nlay_i + nlay_s, nlay_s + 2, -1 518 DO ji = 1, n idx518 DO ji = 1, npti 519 519 jk = jm - nlay_s - 1 520 520 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) … … 522 522 END DO 523 523 524 DO ji = 1, n idx524 DO ji = 1, npti 525 525 ! snow temperatures 526 526 IF( h_s_1d(ji) > 0._wp ) THEN … … 542 542 ! zdti_max is a measure of error, it has to be under zdti_bnd 543 543 zdti_max = 0._wp 544 DO ji = 1, n idx544 DO ji = 1, npti 545 545 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp ) 546 546 zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) ) … … 548 548 549 549 DO jk = 1, nlay_s 550 DO ji = 1, n idx550 DO ji = 1, npti 551 551 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 552 552 zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) … … 555 555 556 556 DO jk = 1, nlay_i 557 DO ji = 1, n idx557 DO ji = 1, npti 558 558 ztmelt_i = -tmut * sz_i_1d(ji,jk) + rt0 559 559 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp ) … … 576 576 ! 10) Fluxes at the interfaces 577 577 !----------------------------- 578 DO ji = 1, n idx578 DO ji = 1, npti 579 579 ! ! surface ice conduction flux 580 580 fc_su(ji) = - isnow(ji) * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji)) & … … 589 589 ! --- diagnose the change in non-solar flux due to surface temperature change --- ! 590 590 IF ( ln_dqns_i ) THEN 591 DO ji = 1, n idx591 DO ji = 1, npti 592 592 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji) 593 593 END DO … … 597 597 ! hfx_dif = Heat flux used to warm/cool ice in W.m-2 598 598 ! zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 599 DO ji = 1, n idx599 DO ji = 1, npti 600 600 zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i + & 601 601 & SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) … … 614 614 615 615 ! --- Diagnostics SIMIP --- ! 616 DO ji = 1, n idx616 DO ji = 1, npti 617 617 !--- Conduction fluxes (positive downwards) 618 618 diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji) … … 645 645 ! 646 646 DO jk = 1, nlay_i ! Sea ice energy of melting 647 DO ji = 1, n idx647 DO ji = 1, npti 648 648 ztmelts = - tmut * sz_i_1d(ji,jk) 649 649 t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point … … 655 655 END DO 656 656 DO jk = 1, nlay_s ! Snow energy of melting 657 DO ji = 1, n idx657 DO ji = 1, npti 658 658 e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) 659 659 END DO
Note: See TracChangeset
for help on using the changeset viewer.
