Changeset 921 for trunk/NEMO/LIM_SRC_3/limthd.F90
- Timestamp:
- 2008-05-13T10:28:52+02:00 (16 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/LIM_SRC_3/limthd.F90 (modified) (43 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/LIM_SRC_3/limthd.F90
r888 r921 39 39 PUBLIC lim_thd ! called by lim_step 40 40 41 !! * Module variables42 43 44 45 46 47 48 41 !! * Module variables 42 REAL(wp) :: & ! constant values 43 epsi20 = 1e-20 , & 44 epsi16 = 1e-16 , & 45 epsi06 = 1e-06 , & 46 epsi04 = 1e-04 , & 47 zzero = 0.e0 , & 48 zone = 1.e0 49 49 50 50 !! * Substitutions … … 59 59 CONTAINS 60 60 61 SUBROUTINE lim_thd 61 SUBROUTINE lim_thd( kt ) 62 62 !!------------------------------------------------------------------- 63 63 !! *** ROUTINE lim_thd *** … … 84 84 !! salinity variations 85 85 !!--------------------------------------------------------------------- 86 INTEGER, INTENT(in) :: kt ! number of iteration 86 87 !! * Local variables 87 88 INTEGER :: ji, jj, jk, jl, nbpb ! nb of icy pts for thermo. cal. … … 90 91 zfric_umin = 5e-03 , & ! lower bound for the friction velocity 91 92 zfric_umax = 2e-02 ! upper bound for the friction velocity 92 93 93 94 REAL(wp) :: & 94 95 zinda , & ! switch for test. the val. of concen. … … 103 104 REAL(wp) :: & 104 105 zareamin 105 106 106 107 REAL(wp), DIMENSION(jpi,jpj) :: & 107 108 zhicifp , & ! ice thickness for outputs … … 112 113 IF( numit == nstart ) CALL lim_thd_init ! Initialization (first time-step only) 113 114 114 WRITE(numout,*) 'limthd : Ice Thermodynamics' 115 WRITE(numout,*) '~~~~~~' 115 IF( kt == nit000 .AND. lwp ) THEN 116 WRITE(numout,*) 'limthd : Ice Thermodynamics' 117 WRITE(numout,*) '~~~~~~' 118 ENDIF 116 119 117 120 IF( numit == nstart ) CALL lim_thd_sal_init ! Initialization (first time-step only) 118 !------------------------------------------------------------------------------!119 ! 1) Initialization of diagnostic variables !120 !------------------------------------------------------------------------------!121 !------------------------------------------------------------------------------! 122 ! 1) Initialization of diagnostic variables ! 123 !------------------------------------------------------------------------------! 121 124 zeps = 1.0e-10 122 125 tatm_ice(:,:) = tatm_ice(:,:) + 273.15 ! convert C to K … … 129 132 130 133 DO jl = 1, jpl 131 DO jk = 1, nlay_i 132 DO jj = 1, jpj 133 DO ji = 1, jpi 134 !Energy of melting q(S,T) [J.m-3] 135 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / area(ji,jj) / & 136 MAX( v_i(ji,jj,jl) , epsi06 ) * nlay_i 137 !0 if no ice and 1 if yes 138 zindb = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_i(ji,jj,jl) ) ) 139 !convert units ! very important that this line is here 140 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb 134 DO jk = 1, nlay_i 135 DO jj = 1, jpj 136 DO ji = 1, jpi 137 !Energy of melting q(S,T) [J.m-3] 138 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / area(ji,jj) / & 139 MAX( v_i(ji,jj,jl) , epsi06 ) * nlay_i 140 !0 if no ice and 1 if yes 141 zindb = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_i(ji,jj,jl) ) ) 142 !convert units ! very important that this line is here 143 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb 144 END DO 141 145 END DO 142 END DO 143 END DO 146 END DO 144 147 END DO 145 148 146 149 DO jl = 1, jpl 147 DO jk = 1, nlay_s 148 DO jj = 1, jpj 149 DO ji = 1, jpi 150 !Energy of melting q(S,T) [J.m-3] 151 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / area(ji,jj) / & 152 MAX( v_s(ji,jj,jl) , epsi06 ) * nlay_s 153 !0 if no ice and 1 if yes 154 zindb = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s(ji,jj,jl) ) ) 155 !convert units ! very important that this line is here 156 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * unit_fac * zindb 150 DO jk = 1, nlay_s 151 DO jj = 1, jpj 152 DO ji = 1, jpi 153 !Energy of melting q(S,T) [J.m-3] 154 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / area(ji,jj) / & 155 MAX( v_s(ji,jj,jl) , epsi06 ) * nlay_s 156 !0 if no ice and 1 if yes 157 zindb = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s(ji,jj,jl) ) ) 158 !convert units ! very important that this line is here 159 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * unit_fac * zindb 160 END DO 157 161 END DO 158 END DO 159 END DO 162 END DO 160 163 END DO 161 164 … … 187 190 fatm(:,:) = 0.e0 188 191 189 ! 2) Partial computation of forcing for the thermodynamic sea ice model. !190 !-----------------------------------------------------------------------------!191 192 ! !CDIR NOVERRCHK193 194 ! !CDIR NOVERRCHK195 192 ! 2) Partial computation of forcing for the thermodynamic sea ice model. ! 193 !-----------------------------------------------------------------------------! 194 195 !CDIR NOVERRCHK 196 DO jj = 1, jpj 197 !CDIR NOVERRCHK 198 DO ji = 1, jpi 196 199 zthsnice = SUM( ht_s(ji,jj,1:jpl) ) + SUM( ht_i(ji,jj,1:jpl) ) 197 200 zindb = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - zthsnice ) ) ) … … 199 202 pfrld(ji,jj) = 1.0 - at_i(ji,jj) 200 203 zinda = 1.0 - MAX( zzero , SIGN( zone , - ( 1.0 - pfrld(ji,jj) ) ) ) 201 202 ! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget203 ! ! practically no "direct lateral ablation"204 !205 ! ! net downward heat flux from the ice to the ocean, expressed as a function of ocean206 ! ! temperature and turbulent mixing (McPhee, 1992)204 205 ! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget 206 ! ! practically no "direct lateral ablation" 207 ! 208 ! ! net downward heat flux from the ice to the ocean, expressed as a function of ocean 209 ! ! temperature and turbulent mixing (McPhee, 1992) 207 210 ! friction velocity 208 211 zfric_u = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin ) … … 211 214 fdtcn(ji,jj) = zindb * rau0 * rcp * 0.006 * zfric_u * ( (sst_m(ji,jj) + rt0) - t_bo(ji,jj) ) 212 215 ! also category dependent 213 ! !-- Energy from the turbulent oceanic heat flux heat flux coming in the lead216 ! !-- Energy from the turbulent oceanic heat flux heat flux coming in the lead 214 217 qdtcn(ji,jj) = zindb * fdtcn(ji,jj) * (1.0 - at_i(ji,jj)) * rdt_ice 215 !216 217 ! still need to be updated : fdtcn !!!!218 ! !-- Lead heat budget (part 1, next one is in limthd_dh219 ! !-- qldif -- (or qldif_1d in 1d routines)218 ! 219 220 ! still need to be updated : fdtcn !!!! 221 ! !-- Lead heat budget (part 1, next one is in limthd_dh 222 ! !-- qldif -- (or qldif_1d in 1d routines) 220 223 zfontn = sprecip(ji,jj) * lfus ! energy of melting 221 224 zfnsol = qns(ji,jj) ! total non solar flux … … 232 235 !false 233 236 zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / & 234 237 MAX( at_i(ji,jj) * rdt_ice , epsi16 ) 235 238 236 239 ! Heat budget of the lead, energy transferred from ice to ocean … … 244 247 ! calculate oceanic heat flux (limthd_dh) 245 248 fbif (ji,jj) = zindb * ( fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi20 ) + fdtcn(ji,jj) ) 246 249 247 250 ! computation of the daily thermodynamic ice production (only needed for output) 248 251 zhicifp(ji,jj) = ht_i(ji,jj,1) * at_i(ji,jj) … … 251 254 END DO 252 255 253 !------------------------------------------------------------------------------!254 ! 3) Select icy points and fulfill arrays for the vectorial grid.255 !------------------------------------------------------------------------------!256 !------------------------------------------------------------------------------! 257 ! 3) Select icy points and fulfill arrays for the vectorial grid. 258 !------------------------------------------------------------------------------! 256 259 257 260 DO jl = 1, jpl !loop over ice categories 258 261 259 WRITE(numout,*) ' lim_thd : transfer to 1D vectors. Category no : ', jl 260 WRITE(numout,*) ' ~~~~~~~~' 262 IF( kt == nit000 .AND. lwp ) THEN 263 WRITE(numout,*) ' lim_thd : transfer to 1D vectors. Category no : ', jl 264 WRITE(numout,*) ' ~~~~~~~~' 265 ENDIF 261 266 262 267 zareamin = 1.0e-10 … … 270 275 ! debug point to follow 271 276 IF ( (ji.eq.jiindx).AND.(jj.eq.jjindx) ) THEN 272 277 jiindex_1d = nbpb 273 278 ENDIF 274 279 END DO 275 280 END DO 276 281 277 !------------------------------------------------------------------------------!278 ! 4) Thermodynamic computation279 !------------------------------------------------------------------------------!282 !------------------------------------------------------------------------------! 283 ! 4) Thermodynamic computation 284 !------------------------------------------------------------------------------! 280 285 281 286 IF( lk_mpp ) CALL mpp_ini_ice(nbpb) … … 283 288 IF (nbpb > 0) THEN ! If there is no ice, do nothing. 284 289 285 !-------------------------286 ! 4.1 Move to 1D arrays287 !-------------------------290 !------------------------- 291 ! 4.1 Move to 1D arrays 292 !------------------------- 288 293 289 294 CALL tab_2d_1d( nbpb, at_i_b (1:nbpb) , at_i , jpi, jpj, npb(1:nbpb) ) … … 330 335 CALL tab_2d_1d( nbpb, qfvbq_1d (1:nbpb) , qfvbq , jpi, jpj, npb(1:nbpb) ) 331 336 332 !--------------------------------333 ! 4.3) Thermodynamic processes334 !--------------------------------335 337 !-------------------------------- 338 ! 4.3) Thermodynamic processes 339 !-------------------------------- 340 336 341 IF ( con_i ) CALL lim_thd_enmelt(1,nbpb) ! computes sea ice energy of melting 337 342 IF ( con_i ) CALL lim_thd_glohec( qt_i_in , qt_s_in , & 338 339 340 343 q_i_layer_in , 1 , nbpb , jl ) 344 345 !---------------------------------! 341 346 CALL lim_thd_dif(1,nbpb,jl) ! Ice/Snow Temperature profile ! 342 347 !---------------------------------! 343 348 344 349 CALL lim_thd_enmelt(1,nbpb) ! computes sea ice energy of melting 345 350 ! compulsory for limthd_dh 346 351 347 352 IF ( con_i ) CALL lim_thd_glohec( qt_i_fin , qt_s_fin , & 348 353 q_i_layer_fin , 1 , nbpb , jl ) 349 354 IF ( con_i ) CALL lim_thd_con_dif( 1 , nbpb , jl ) 350 355 351 356 !---------------------------------! 352 357 CALL lim_thd_dh(1,nbpb,jl) ! Ice/Snow thickness ! 353 354 355 358 !---------------------------------! 359 360 !---------------------------------! 356 361 CALL lim_thd_ent(1,nbpb,jl) ! Ice/Snow enthalpy remapping ! 357 358 359 362 !---------------------------------! 363 364 !---------------------------------! 360 365 CALL lim_thd_sal(1,nbpb) ! Ice salinity computation ! 361 362 363 ! CALL lim_thd_enmelt(1,nbpb) ! computes sea ice energy of melting366 !---------------------------------! 367 368 ! CALL lim_thd_enmelt(1,nbpb) ! computes sea ice energy of melting 364 369 IF ( con_i ) CALL lim_thd_glohec( qt_i_fin, qt_s_fin, & 365 370 q_i_layer_fin , 1 , nbpb , jl ) 366 371 IF ( con_i ) CALL lim_thd_con_dh ( 1 , nbpb , jl ) 367 372 368 !--------------------------------369 ! 4.4) Move 1D to 2D vectors370 !--------------------------------373 !-------------------------------- 374 ! 4.4) Move 1D to 2D vectors 375 !-------------------------------- 371 376 372 377 CALL tab_1d_2d( nbpb, at_i , npb, at_i_b (1:nbpb), jpi, jpj ) … … 416 421 !+++++ 417 422 418 IF( lk_mpp ) CALL mpp_comm_free(ncomm_ice) !RB necessary ??423 IF( lk_mpp ) CALL mpp_comm_free(ncomm_ice) !RB necessary ?? 419 424 ENDIF ! nbpb 420 425 421 426 END DO ! jl 422 427 423 !------------------------------------------------------------------------------!424 ! 5) Global variables, diagnostics425 !------------------------------------------------------------------------------!428 !------------------------------------------------------------------------------! 429 ! 5) Global variables, diagnostics 430 !------------------------------------------------------------------------------! 426 431 427 432 !------------------------ … … 431 436 ! Enthalpies are global variables we have to readjust the units 432 437 DO jl = 1, jpl 433 DO jk = 1, nlay_i434 DO jj = 1, jpj435 DO ji = 1, jpi436 ! Change dimensions437 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac438 439 ! Multiply by volume, divide by nlayers so that heat content in 10^9 Joules440 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * &441 442 443 END DO !ji444 END DO !jj445 END DO !jk438 DO jk = 1, nlay_i 439 DO jj = 1, jpj 440 DO ji = 1, jpi 441 ! Change dimensions 442 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac 443 444 ! Multiply by volume, divide by nlayers so that heat content in 10^9 Joules 445 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * & 446 area(ji,jj) * a_i(ji,jj,jl) * & 447 ht_i(ji,jj,jl) / nlay_i 448 END DO !ji 449 END DO !jj 450 END DO !jk 446 451 END DO !jl 447 452 … … 459 464 ! Multiply by volume, so that heat content in 10^9 Joules 460 465 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * & 461 466 a_i(ji,jj,jl) * ht_s(ji,jj,jl) / nlay_s 462 467 END DO !ji 463 468 END DO !jj … … 513 518 END SUBROUTINE lim_thd 514 519 515 !===============================================================================520 !=============================================================================== 516 521 517 522 SUBROUTINE lim_thd_glohec(eti,ets,etilayer,kideb,kiut,jl) … … 552 557 DO ji = kideb, kiut 553 558 etilayer(ji,jk) = q_i_b(ji,jk) & 554 559 * ht_i_b(ji) / nlay_i 555 560 eti(ji,jl) = eti(ji,jl) + etilayer(ji,jk) 556 561 END DO … … 567 572 WRITE(numout,*) ' qt_s_in : ', ets(jiindex_1d,jl) / rdt_ice 568 573 WRITE(numout,*) ' qt_in : ', ( eti(jiindex_1d,jl) + & 569 574 ets(jiindex_1d,jl) ) / rdt_ice 570 575 571 576 END SUBROUTINE lim_thd_glohec 572 577 573 !===============================================================================578 !=============================================================================== 574 579 575 580 SUBROUTINE lim_thd_con_dif(kideb,kiut,jl) … … 594 599 INTEGER :: & 595 600 numce !: number of points for which conservation 596 601 ! is violated 597 602 INTEGER :: & 598 603 ji,jk, & !: loop indices … … 602 607 max_cons_err = 1.0 603 608 max_surf_err = 0.001 604 609 605 610 !-------------------------- 606 611 ! Increment of energy … … 608 613 ! global 609 614 DO ji = kideb, kiut 610 611 615 dq_i(ji,jl) = qt_i_fin(ji,jl) - qt_i_in(ji,jl) & 616 + qt_s_fin(ji,jl) - qt_s_in(ji,jl) 612 617 END DO 613 618 ! layer by layer … … 619 624 620 625 DO ji = kideb, kiut 621 622 623 624 625 qnsr_ice_1d(ji) + & ! atm non solar626 (1.0-i0(ji))*qsr_ice_1d(ji) ! atm solar627 628 629 626 zji = MOD( npb(ji) - 1, jpi ) + 1 627 zjj = ( npb(ji) - 1 ) / jpi + 1 628 629 fatm(ji,jl) = & 630 qnsr_ice_1d(ji) + & ! atm non solar 631 (1.0-i0(ji))*qsr_ice_1d(ji) ! atm solar 632 633 sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji) + qsr_ice_1d(ji)*i0(ji) & 634 - fstroc(zji,zjj,jl) 630 635 END DO 631 636 … … 635 640 636 641 DO ji = kideb, kiut 637 642 cons_error(ji,jl) = ABS( dq_i(ji,jl) / rdt_ice + sum_fluxq(ji,jl) ) 638 643 END DO 639 644 … … 641 646 meance = 0.0 642 647 DO ji = kideb, kiut 643 644 645 646 648 IF ( cons_error(ji,jl) .GT. max_cons_err ) THEN 649 numce = numce + 1 650 meance = meance + cons_error(ji,jl) 651 ENDIF 647 652 ENDDO 648 653 IF (numce .GT. 0 ) meance = meance / numce … … 651 656 WRITE(numout,*) ' After lim_thd_dif, category : ', jl 652 657 WRITE(numout,*) ' Mean conservation error on big error points ', meance, & 653 numit658 numit 654 659 WRITE(numout,*) ' Number of points where there is a cons err gt than c.e. : ', numce, numit 655 660 … … 663 668 surf_error(ji,jl) = ABS ( fatm(ji,jl) - fc_su(ji) ) 664 669 IF ( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. & 665 670 max_surf_err ) ) THEN 666 671 numce = numce + 1 667 672 meance = meance + surf_error(ji,jl) … … 685 690 DO ji = kideb, kiut 686 691 IF ( ( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. max_surf_err ) ) .OR. & 687 688 zji = MOD( npb(ji) - 1, jpi ) + 1689 zjj = ( npb(ji) - 1 ) / jpi + 1690 691 WRITE(numout,*) ' alerte 1 '692 WRITE(numout,*) ' Untolerated conservation / surface error after '693 WRITE(numout,*) ' heat diffusion in the ice '694 WRITE(numout,*) ' Category : ', jl695 WRITE(numout,*) ' zji , zjj : ', zji, zjj696 WRITE(numout,*) ' lat, lon : ', gphit(zji,zjj), glamt(zji,zjj)697 WRITE(numout,*) ' cons_error : ', cons_error(ji,jl)698 WRITE(numout,*) ' surf_error : ', surf_error(ji,jl)699 WRITE(numout,*) ' dq_i : ', - dq_i(ji,jl) / rdt_ice700 WRITE(numout,*) ' Fdt : ', sum_fluxq(ji,jl)701 WRITE(numout,*)702 ! WRITE(numout,*) ' qt_i_in : ', qt_i_in(ji,jl)703 ! WRITE(numout,*) ' qt_s_in : ', qt_s_in(ji,jl)704 ! WRITE(numout,*) ' qt_i_fin : ', qt_i_fin(ji,jl)705 ! WRITE(numout,*) ' qt_s_fin : ', qt_s_fin(ji,jl)706 ! WRITE(numout,*) ' qt : ', qt_i_fin(ji,jl) + &707 ! qt_s_fin(ji,jl)708 WRITE(numout,*) ' ht_i : ', ht_i_b(ji)709 WRITE(numout,*) ' ht_s : ', ht_s_b(ji)710 WRITE(numout,*) ' t_su : ', t_su_b(ji)711 WRITE(numout,*) ' t_s : ', t_s_b(ji,1)712 WRITE(numout,*) ' t_i : ', t_i_b(ji,1:nlay_i)713 WRITE(numout,*) ' t_bo : ', t_bo_b(ji)714 WRITE(numout,*) ' q_i : ', q_i_b(ji,1:nlay_i)715 WRITE(numout,*) ' s_i : ', s_i_b(ji,1:nlay_i)716 WRITE(numout,*) ' tmelts : ', rtt - tmut*s_i_b(ji,1:nlay_i)717 WRITE(numout,*)718 WRITE(numout,*) ' Fluxes '719 WRITE(numout,*) ' ~~~~~~ '720 WRITE(numout,*) ' fatm : ', fatm(ji,jl)721 WRITE(numout,*) ' fc_su : ', fc_su (ji)722 WRITE(numout,*) ' fstr_inice : ', qsr_ice_1d(ji)*i0(ji)723 WRITE(numout,*) ' fc_bo : ', - fc_bo_i (ji)724 WRITE(numout,*) ' foc : ', fbif_1d(ji)725 WRITE(numout,*) ' fstroc : ', fstroc (zji,zjj,jl)726 WRITE(numout,*) ' i0 : ', i0(ji)727 WRITE(numout,*) ' qsr_ice : ', (1.0-i0(ji))*qsr_ice_1d(ji)728 WRITE(numout,*) ' qns_ice : ', qnsr_ice_1d(ji)729 WRITE(numout,*) ' Conduction fluxes : '730 WRITE(numout,*) ' fc_s : ', fc_s(ji,0:nlay_s)731 WRITE(numout,*) ' fc_i : ', fc_i(ji,0:nlay_i)732 WRITE(numout,*)733 WRITE(numout,*) ' Layer by layer ... '734 WRITE(numout,*) ' dq_snow : ', ( qt_s_fin(ji,jl) - &735 736 737 WRITE(numout,*) ' dfc_snow : ', fc_s(ji,1) - &738 739 DO jk = 1, nlay_i740 WRITE(numout,*) ' layer : ', jk741 WRITE(numout,*) ' dq_ice : ', dq_i_layer(ji,jk) / rdt_ice742 WRITE(numout,*) ' radab : ', radab(ji,jk)743 WRITE(numout,*) ' dfc_i : ', fc_i(ji,jk) - &744 745 WRITE(numout,*) ' tot f : ', fc_i(ji,jk) - &746 747 END DO692 ( cons_error(ji,jl) .GT. max_cons_err ) ) THEN 693 zji = MOD( npb(ji) - 1, jpi ) + 1 694 zjj = ( npb(ji) - 1 ) / jpi + 1 695 696 WRITE(numout,*) ' alerte 1 ' 697 WRITE(numout,*) ' Untolerated conservation / surface error after ' 698 WRITE(numout,*) ' heat diffusion in the ice ' 699 WRITE(numout,*) ' Category : ', jl 700 WRITE(numout,*) ' zji , zjj : ', zji, zjj 701 WRITE(numout,*) ' lat, lon : ', gphit(zji,zjj), glamt(zji,zjj) 702 WRITE(numout,*) ' cons_error : ', cons_error(ji,jl) 703 WRITE(numout,*) ' surf_error : ', surf_error(ji,jl) 704 WRITE(numout,*) ' dq_i : ', - dq_i(ji,jl) / rdt_ice 705 WRITE(numout,*) ' Fdt : ', sum_fluxq(ji,jl) 706 WRITE(numout,*) 707 ! WRITE(numout,*) ' qt_i_in : ', qt_i_in(ji,jl) 708 ! WRITE(numout,*) ' qt_s_in : ', qt_s_in(ji,jl) 709 ! WRITE(numout,*) ' qt_i_fin : ', qt_i_fin(ji,jl) 710 ! WRITE(numout,*) ' qt_s_fin : ', qt_s_fin(ji,jl) 711 ! WRITE(numout,*) ' qt : ', qt_i_fin(ji,jl) + & 712 ! qt_s_fin(ji,jl) 713 WRITE(numout,*) ' ht_i : ', ht_i_b(ji) 714 WRITE(numout,*) ' ht_s : ', ht_s_b(ji) 715 WRITE(numout,*) ' t_su : ', t_su_b(ji) 716 WRITE(numout,*) ' t_s : ', t_s_b(ji,1) 717 WRITE(numout,*) ' t_i : ', t_i_b(ji,1:nlay_i) 718 WRITE(numout,*) ' t_bo : ', t_bo_b(ji) 719 WRITE(numout,*) ' q_i : ', q_i_b(ji,1:nlay_i) 720 WRITE(numout,*) ' s_i : ', s_i_b(ji,1:nlay_i) 721 WRITE(numout,*) ' tmelts : ', rtt - tmut*s_i_b(ji,1:nlay_i) 722 WRITE(numout,*) 723 WRITE(numout,*) ' Fluxes ' 724 WRITE(numout,*) ' ~~~~~~ ' 725 WRITE(numout,*) ' fatm : ', fatm(ji,jl) 726 WRITE(numout,*) ' fc_su : ', fc_su (ji) 727 WRITE(numout,*) ' fstr_inice : ', qsr_ice_1d(ji)*i0(ji) 728 WRITE(numout,*) ' fc_bo : ', - fc_bo_i (ji) 729 WRITE(numout,*) ' foc : ', fbif_1d(ji) 730 WRITE(numout,*) ' fstroc : ', fstroc (zji,zjj,jl) 731 WRITE(numout,*) ' i0 : ', i0(ji) 732 WRITE(numout,*) ' qsr_ice : ', (1.0-i0(ji))*qsr_ice_1d(ji) 733 WRITE(numout,*) ' qns_ice : ', qnsr_ice_1d(ji) 734 WRITE(numout,*) ' Conduction fluxes : ' 735 WRITE(numout,*) ' fc_s : ', fc_s(ji,0:nlay_s) 736 WRITE(numout,*) ' fc_i : ', fc_i(ji,0:nlay_i) 737 WRITE(numout,*) 738 WRITE(numout,*) ' Layer by layer ... ' 739 WRITE(numout,*) ' dq_snow : ', ( qt_s_fin(ji,jl) - & 740 qt_s_in(ji,jl) ) & 741 / rdt_ice 742 WRITE(numout,*) ' dfc_snow : ', fc_s(ji,1) - & 743 fc_s(ji,0) 744 DO jk = 1, nlay_i 745 WRITE(numout,*) ' layer : ', jk 746 WRITE(numout,*) ' dq_ice : ', dq_i_layer(ji,jk) / rdt_ice 747 WRITE(numout,*) ' radab : ', radab(ji,jk) 748 WRITE(numout,*) ' dfc_i : ', fc_i(ji,jk) - & 749 fc_i(ji,jk-1) 750 WRITE(numout,*) ' tot f : ', fc_i(ji,jk) - & 751 fc_i(ji,jk-1) - radab(ji,jk) 752 END DO 748 753 749 754 ENDIF 750 755 751 756 END DO 752 757 753 758 END SUBROUTINE lim_thd_con_dif 754 759 755 !==============================================================================760 !============================================================================== 756 761 757 762 SUBROUTINE lim_thd_con_dh(kideb,kiut,jl) … … 775 780 INTEGER :: & 776 781 numce !: number of points for which conservation 777 782 ! is violated 778 783 INTEGER :: ji, zji, zjj ! loop indices 779 784 !!--------------------------------------------------------------------- 780 785 781 786 max_cons_err = 1.0 782 787 783 788 !-------------------------- 784 789 ! Increment of energy … … 786 791 ! global 787 792 DO ji = kideb, kiut 788 789 793 dq_i(ji,jl) = qt_i_fin(ji,jl) - qt_i_in(ji,jl) & 794 + qt_s_fin(ji,jl) - qt_s_in(ji,jl) 790 795 END DO 791 796 ! layer by layer … … 797 802 798 803 DO ji = kideb, kiut 799 800 801 802 803 qnsr_ice_1d(ji) + & ! atm non solar804 ! (1.0-i0(ji))*qsr_ice_1d(ji) ! atm solar805 qsr_ice_1d(ji) ! atm solar806 807 808 809 804 zji = MOD( npb(ji) - 1, jpi ) + 1 805 zjj = ( npb(ji) - 1 ) / jpi + 1 806 807 fatm(ji,jl) = & 808 qnsr_ice_1d(ji) + & ! atm non solar 809 ! (1.0-i0(ji))*qsr_ice_1d(ji) ! atm solar 810 qsr_ice_1d(ji) ! atm solar 811 812 sum_fluxq(ji,jl) = fatm(ji,jl) + fbif_1d(ji) - ftotal_fin(ji) & 813 - fstroc(zji,zjj,jl) 814 cons_error(ji,jl) = ABS( dq_i(ji,jl) / rdt_ice + sum_fluxq(ji,jl) ) 810 815 END DO 811 816 … … 815 820 816 821 DO ji = kideb, kiut 817 822 cons_error(ji,jl) = ABS( dq_i(ji,jl) / rdt_ice + sum_fluxq(ji,jl) ) 818 823 END DO 819 824 … … 821 826 meance = 0.0 822 827 DO ji = kideb, kiut 823 824 825 826 828 IF ( cons_error(ji,jl) .GT. max_cons_err ) THEN 829 numce = numce + 1 830 meance = meance + cons_error(ji,jl) 831 ENDIF 827 832 ENDDO 828 833 IF (numce .GT. 0 ) meance = meance / numce … … 833 838 WRITE(numout,*) ' After lim_thd_ent, category : ', jl 834 839 WRITE(numout,*) ' Mean conservation error on big error points ', meance, & 835 numit840 numit 836 841 WRITE(numout,*) ' Number of points where there is a cons err gt than 0.1 W/m2 : ', numce, numit 837 842 … … 842 847 DO ji = kideb, kiut 843 848 IF ( cons_error(ji,jl) .GT. max_cons_err ) THEN 844 zji = MOD( npb(ji) - 1, jpi ) + 1845 zjj = ( npb(ji) - 1 ) / jpi + 1846 847 WRITE(numout,*) ' alerte 1 - category : ', jl848 WRITE(numout,*) ' Untolerated conservation error after limthd_ent '849 WRITE(numout,*) ' zji , zjj : ', zji, zjj850 WRITE(numout,*) ' lat, lon : ', gphit(zji,zjj), glamt(zji,zjj)851 WRITE(numout,*) ' * '852 WRITE(numout,*) ' Ftotal : ', sum_fluxq(ji,jl)853 WRITE(numout,*) ' dq_t : ', - dq_i(ji,jl) / rdt_ice854 WRITE(numout,*) ' dq_i : ', - ( qt_i_fin(ji,jl) - qt_i_in(ji,jl) ) / rdt_ice855 WRITE(numout,*) ' dq_s : ', - ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) / rdt_ice856 WRITE(numout,*) ' cons_error : ', cons_error(ji,jl)857 WRITE(numout,*) ' * '858 WRITE(numout,*) ' Fluxes --- : '859 WRITE(numout,*) ' fatm : ', fatm(ji,jl)860 WRITE(numout,*) ' foce : ', fbif_1d(ji)861 WRITE(numout,*) ' fres : ', ftotal_fin(ji)862 WRITE(numout,*) ' fhbri : ', fhbricat(zji,zjj,jl)863 WRITE(numout,*) ' * '864 WRITE(numout,*) ' Heat contents --- : '865 WRITE(numout,*) ' qt_s_in : ', qt_s_in(ji,jl) / rdt_ice866 WRITE(numout,*) ' qt_i_in : ', qt_i_in(ji,jl) / rdt_ice867 WRITE(numout,*) ' qt_in : ', ( qt_i_in(ji,jl) + &868 869 WRITE(numout,*) ' qt_s_fin : ', qt_s_fin(ji,jl) / rdt_ice870 WRITE(numout,*) ' qt_i_fin : ', qt_i_fin(ji,jl) / rdt_ice871 WRITE(numout,*) ' qt_fin : ', ( qt_i_fin(ji,jl) + &872 873 WRITE(numout,*) ' * '874 WRITE(numout,*) ' Ice variables --- : '875 WRITE(numout,*) ' ht_i : ', ht_i_b(ji)876 WRITE(numout,*) ' ht_s : ', ht_s_b(ji)877 WRITE(numout,*) ' dh_s_tot : ', dh_s_tot(ji)878 WRITE(numout,*) ' dh_snowice: ', dh_snowice(ji)879 WRITE(numout,*) ' dh_i_surf : ', dh_i_surf(ji)880 WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)849 zji = MOD( npb(ji) - 1, jpi ) + 1 850 zjj = ( npb(ji) - 1 ) / jpi + 1 851 852 WRITE(numout,*) ' alerte 1 - category : ', jl 853 WRITE(numout,*) ' Untolerated conservation error after limthd_ent ' 854 WRITE(numout,*) ' zji , zjj : ', zji, zjj 855 WRITE(numout,*) ' lat, lon : ', gphit(zji,zjj), glamt(zji,zjj) 856 WRITE(numout,*) ' * ' 857 WRITE(numout,*) ' Ftotal : ', sum_fluxq(ji,jl) 858 WRITE(numout,*) ' dq_t : ', - dq_i(ji,jl) / rdt_ice 859 WRITE(numout,*) ' dq_i : ', - ( qt_i_fin(ji,jl) - qt_i_in(ji,jl) ) / rdt_ice 860 WRITE(numout,*) ' dq_s : ', - ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) / rdt_ice 861 WRITE(numout,*) ' cons_error : ', cons_error(ji,jl) 862 WRITE(numout,*) ' * ' 863 WRITE(numout,*) ' Fluxes --- : ' 864 WRITE(numout,*) ' fatm : ', fatm(ji,jl) 865 WRITE(numout,*) ' foce : ', fbif_1d(ji) 866 WRITE(numout,*) ' fres : ', ftotal_fin(ji) 867 WRITE(numout,*) ' fhbri : ', fhbricat(zji,zjj,jl) 868 WRITE(numout,*) ' * ' 869 WRITE(numout,*) ' Heat contents --- : ' 870 WRITE(numout,*) ' qt_s_in : ', qt_s_in(ji,jl) / rdt_ice 871 WRITE(numout,*) ' qt_i_in : ', qt_i_in(ji,jl) / rdt_ice 872 WRITE(numout,*) ' qt_in : ', ( qt_i_in(ji,jl) + & 873 qt_s_in(ji,jl) ) / rdt_ice 874 WRITE(numout,*) ' qt_s_fin : ', qt_s_fin(ji,jl) / rdt_ice 875 WRITE(numout,*) ' qt_i_fin : ', qt_i_fin(ji,jl) / rdt_ice 876 WRITE(numout,*) ' qt_fin : ', ( qt_i_fin(ji,jl) + & 877 qt_s_fin(ji,jl) ) / rdt_ice 878 WRITE(numout,*) ' * ' 879 WRITE(numout,*) ' Ice variables --- : ' 880 WRITE(numout,*) ' ht_i : ', ht_i_b(ji) 881 WRITE(numout,*) ' ht_s : ', ht_s_b(ji) 882 WRITE(numout,*) ' dh_s_tot : ', dh_s_tot(ji) 883 WRITE(numout,*) ' dh_snowice: ', dh_snowice(ji) 884 WRITE(numout,*) ' dh_i_surf : ', dh_i_surf(ji) 885 WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji) 881 886 882 887 ENDIF 883 888 884 889 END DO 885 890 886 891 END SUBROUTINE lim_thd_con_dh 887 !==============================================================================892 !============================================================================== 888 893 889 894 SUBROUTINE lim_thd_enmelt(kideb,kiut) … … 899 904 INTEGER, INTENT(in) :: & 900 905 kideb, kiut !: bounds for the spatial loop 901 906 902 907 REAL(wp) :: & !: goes to trash 903 908 ztmelts , & !: sea ice freezing point in K … … 916 921 ztmelts = - tmut * s_i_b(ji,jk) + rtt 917 922 q_i_b(ji,jk) = rhoic*( cpic * ( ztmelts - t_i_b(ji,jk) ) & 918 919 923 + lfus * ( 1.0 - (ztmelts-rtt)/MIN((t_i_b(ji,jk)-rtt),-zeps) ) & 924 - rcp * ( ztmelts-rtt ) ) 920 925 END DO !ji 921 926 END DO !jk … … 930 935 END SUBROUTINE lim_thd_enmelt 931 936 932 !==============================================================================937 !============================================================================== 933 938 934 939 SUBROUTINE lim_thd_init … … 954 959 & kappa_i, nconv_i_thd, maxer_i_thd, thcon_i_swi 955 960 !!------------------------------------------------------------------- 956 961 957 962 ! Define the initial parameters 958 963 ! ------------------------- … … 990 995 WRITE(numout,*) 991 996 ENDIF 992 997 993 998 rcdsn = hakdif * rcdsn 994 999 rcdic = hakdif * rcdic 995 1000 996 1001 997 1002 END SUBROUTINE lim_thd_init
Note: See TracChangeset
for help on using the changeset viewer.
