- 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.F90
r4147 r4161 11 11 !! 3.3 ! 2010-11 (G. Madec) corrected snow melting heat (due to factor betas) 12 12 !! 4.0 ! 2011-02 (G. Madec) dynamical allocation 13 !! - ! 2012-05 (C. Rousset) add penetration solar flux 13 14 !!---------------------------------------------------------------------- 14 15 #if defined key_lim3 … … 40 41 USE prtctl ! Print control 41 42 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 43 USE timing ! Timing 42 44 43 45 IMPLICIT NONE … … 92 94 REAL(wp) :: zfntlat, zpareff, zareamin, zcoef ! - - 93 95 REAL(wp), POINTER, DIMENSION(:,:) :: zqlbsbq ! link with lead energy budget qldif 96 REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 97 REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 94 98 !!------------------------------------------------------------------- 99 IF( nn_timing == 1 ) CALL timing_start('limthd') 95 100 96 101 CALL wrk_alloc( jpi, jpj, zqlbsbq ) 97 102 103 ! ------------------------------- 104 !- check conservation (C Rousset) 105 IF (ln_limdiahsb) THEN 106 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 107 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 108 zchk_fw_b = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 109 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 110 ENDIF 111 !- check conservation (C Rousset) 112 ! ------------------------------- 113 98 114 !------------------------------------------------------------------------------! 99 115 ! 1) Initialization of diagnostic variables ! … … 109 125 DO ji = 1, jpi 110 126 !Energy of melting q(S,T) [J.m-3] 111 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi06 ) ) * nlay_i127 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi06 ) ) * REAL( nlay_i ) 112 128 !0 if no ice and 1 if yes 113 129 zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - ht_i(ji,jj,jl) ) ) … … 121 137 DO ji = 1, jpi 122 138 !Energy of melting q(S,T) [J.m-3] 123 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi06 ) ) * nlay_s139 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi06 ) ) * REAL( nlay_s ) 124 140 !0 if no ice and 1 if yes 125 141 zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - ht_s(ji,jj,jl) ) ) … … 134 150 ! 1.3) Set some dummies to 0 135 151 !----------------------------- 136 rdvosif(:,:) = 0.e0 ! variation of ice volume at surface137 rdvobif(:,:) = 0.e0 ! variation of ice volume at bottom138 fdvolif(:,:) = 0.e0 ! total variation of ice volume139 rdvonif(:,:) = 0.e0 ! lateral variation of ice volume140 fstric (:,:) = 0.e0 ! part of solar radiation transmitted through the ice141 ffltbif(:,:) = 0.e0 ! linked with fstric142 qfvbq (:,:) = 0.e0 ! linked with fstric143 rdm_snw(:,:) = 0.e0 ! variation of snow mass per unit area144 rdm_ice(:,:) = 0.e0 ! variation of ice mass per unit area145 hicifp (:,:) = 0.e0 ! daily thermodynamic ice production.146 sfx_bri(:,:) = 0.e0 ! brine flux contribution to salt flux to the ocean147 fhbri (:,:) = 0.e0 ! brine flux contribution to heat flux to the ocean148 sfx_thd(:,:) = 0.e0 ! equivalent salt flux to the ocean due to ice/growth decay152 !clem rdvosif(:,:) = 0.e0 ! variation of ice volume at surface 153 !clem rdvobif(:,:) = 0.e0 ! variation of ice volume at bottom 154 !clem fdvolif(:,:) = 0.e0 ! total variation of ice volume 155 !clem rdvonif(:,:) = 0.e0 ! lateral variation of ice volume 156 !clem fstric (:,:) = 0.e0 ! part of solar radiation transmitted through the ice 157 !clem ffltbif(:,:) = 0.e0 ! linked with fstric 158 !clem qfvbq (:,:) = 0.e0 ! linked with fstric 159 !clem rdm_snw(:,:) = 0.e0 ! variation of snow mass per unit area 160 !clem rdm_ice(:,:) = 0.e0 ! variation of ice mass per unit area 161 !clem hicifp (:,:) = 0.e0 ! daily thermodynamic ice production. 162 !clem sfx_bri(:,:) = 0.e0 ! brine flux contribution to salt flux to the ocean 163 !clem fhbri (:,:) = 0.e0 ! brine flux contribution to heat flux to the ocean 164 !clem sfx_thd(:,:) = 0.e0 ! equivalent salt flux to the ocean due to ice/growth decay 149 165 150 166 !----------------------------------- … … 165 181 !CDIR NOVERRCHK 166 182 DO ji = 1, jpi 167 zthsnice = SUM( ht_s(ji,jj,1:jpl) ) + SUM( ht_i(ji,jj,1:jpl) )168 zindb = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - zthsnice ) ) )169 183 phicif(ji,jj) = vt_i(ji,jj) 170 184 pfrld(ji,jj) = 1.0 - at_i(ji,jj) 171 zinda = 1.0 - MAX( zzero , SIGN( zone , - ( 1.0 - pfrld(ji,jj) ) ) )185 zinda = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - at_i(ji,jj) ) ) ) 172 186 ! 173 187 ! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget … … 180 194 181 195 ! here the drag will depend on ice thickness and type (0.006) 182 fdtcn(ji,jj) = zind b* rau0 * rcp * 0.006 * zfric_u * ( (sst_m(ji,jj) + rt0) - t_bo(ji,jj) )196 fdtcn(ji,jj) = zinda * rau0 * rcp * 0.006 * zfric_u * ( (sst_m(ji,jj) + rt0) - t_bo(ji,jj) ) 183 197 ! also category dependent 184 198 ! !-- Energy from the turbulent oceanic heat flux heat flux coming in the lead 185 qdtcn(ji,jj) = zind b* fdtcn(ji,jj) * (1.0 - at_i(ji,jj)) * rdt_ice199 qdtcn(ji,jj) = zinda * fdtcn(ji,jj) * (1.0 - at_i(ji,jj)) * rdt_ice 186 200 ! 187 201 ! !-- Lead heat budget, qldif (part 1, next one is in limthd_dh) 188 202 ! ! caution: exponent betas used as more snow can fallinto leads 189 203 qldif(ji,jj) = tms(ji,jj) * rdt_ice * ( & 190 & pfrld(ji,jj) * ( qsr(ji,jj) & ! solar heat204 & pfrld(ji,jj) * ( qsr(ji,jj) * oatte(ji,jj) & ! solar heat + clem modif 191 205 & + qns(ji,jj) & ! non solar heat 192 206 & + fdtcn(ji,jj) & ! turbulent ice-ocean heat 193 & + fsbbq(ji,jj) * ( 1.0 - zind b) ) & ! residual heat from previous step207 & + fsbbq(ji,jj) * ( 1.0 - zinda ) ) & ! residual heat from previous step 194 208 & - pfrld(ji,jj)**betas * sprecip(ji,jj) * lfus ) ! latent heat of sprecip melting 195 209 ! … … 206 220 ! 207 221 ! Energy needed to bring ocean surface layer until its freezing (qcmif, limflx) 208 qcmif (ji,jj) = rau0 * rcp * fse3t(ji,jj,1) * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) ) * ( 1. - zinda )222 qcmif (ji,jj) = rau0 * rcp * fse3t(ji,jj,1) * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) ) 209 223 ! 210 224 ! oceanic heat flux (limthd_dh) 211 fbif (ji,jj) = zind b* ( fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi20 ) + fdtcn(ji,jj) )225 fbif (ji,jj) = zinda * ( fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi20 ) + fdtcn(ji,jj) ) 212 226 ! 213 227 END DO … … 294 308 CALL tab_2d_1d( nbpb, qfvbq_1d (1:nbpb), qfvbq , jpi, jpj, npb(1:nbpb) ) 295 309 310 CALL tab_2d_1d( nbpb, iatte_1d (1:nbpb), iatte , jpi, jpj, npb(1:nbpb) ) ! clem modif 311 CALL tab_2d_1d( nbpb, oatte_1d (1:nbpb), oatte , jpi, jpj, npb(1:nbpb) ) ! clem modif 296 312 !-------------------------------- 297 313 ! 4.3) Thermodynamic processes … … 411 427 ! 5.4) Diagnostic thermodynamic growth rates 412 428 !-------------------------------------------- 413 d_v_i_thd(:,:,:) = v_i (:,:,:) - old_v_i(:,:,:) ! ice volumes414 429 !clem@useless d_v_i_thd(:,:,:) = v_i (:,:,:) - old_v_i(:,:,:) ! ice volumes 430 !clem@mv-to-itd dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday 415 431 416 432 IF( con_i ) fbif(:,:) = fbif(:,:) + zqlbsbq(:,:) … … 448 464 ENDIF 449 465 ! 466 ! ------------------------------- 467 !- check conservation (C Rousset) 468 IF (ln_limdiahsb) THEN 469 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 470 zchk_fw = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 471 472 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice 473 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 474 475 zchk_vmin = glob_min(v_i) 476 zchk_amax = glob_max(SUM(a_i,dim=3)) 477 zchk_amin = glob_min(a_i) 478 479 IF(lwp) THEN 480 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limthd) = ',(zchk_v_i * rday) 481 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limthd) = ',(zchk_smv * rday) 482 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limthd) = ',(zchk_vmin * 1.e-3) 483 IF ( zchk_amax > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limthd) = ',zchk_amax 484 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limthd) = ',zchk_amin 485 ENDIF 486 ENDIF 487 !- check conservation (C Rousset) 488 ! ------------------------------- 489 ! 450 490 CALL wrk_dealloc( jpi, jpj, zqlbsbq ) 451 491 ! 492 IF( nn_timing == 1 ) CALL timing_stop('limthd') 452 493 END SUBROUTINE lim_thd 453 494 … … 472 513 DO jk = 1, nlay_i ! total q over all layers, ice [J.m-2] 473 514 DO ji = kideb, kiut 474 etilayer(ji,jk) = q_i_b(ji,jk) * ht_i_b(ji) / nlay_i515 etilayer(ji,jk) = q_i_b(ji,jk) * ht_i_b(ji) / REAL( nlay_i ) 475 516 eti (ji,jl) = eti(ji,jl) + etilayer(ji,jk) 476 517 END DO 477 518 END DO 478 519 DO ji = kideb, kiut ! total q over all layers, snow [J.m-2] 479 ets(ji,jl) = ets(ji,jl) + q_s_b(ji,1) * ht_s_b(ji) / nlay_s520 ets(ji,jl) = ets(ji,jl) + q_s_b(ji,1) * ht_s_b(ji) / REAL( nlay_s ) 480 521 END DO 481 522 ! … … 498 539 499 540 INTEGER :: ji, jk ! loop indices 500 INTEGER :: zji, zjj541 INTEGER :: ii, ij 501 542 INTEGER :: numce ! number of points for which conservation is violated 502 543 REAL(wp) :: meance ! mean conservation error … … 521 562 !---------------------------------------- 522 563 DO ji = kideb, kiut 523 zji = MOD( npb(ji) - 1 , jpi ) + 1524 zjj = ( npb(ji) - 1 ) / jpi + 1564 ii = MOD( npb(ji) - 1 , jpi ) + 1 565 ij = ( npb(ji) - 1 ) / jpi + 1 525 566 fatm (ji,jl) = qnsr_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) 526 sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji) + qsr_ice_1d(ji) * i0(ji) - fstroc( zji,zjj,jl)567 sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji) + qsr_ice_1d(ji) * i0(ji) - fstroc(ii,ij,jl) 527 568 END DO 528 569 … … 579 620 IF ( ( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. max_surf_err ) ) .OR. & 580 621 ( cons_error(ji,jl) .GT. max_cons_err ) ) THEN 581 zji = MOD( npb(ji) - 1, jpi ) + 1582 zjj = ( npb(ji) - 1 ) / jpi + 1622 ii = MOD( npb(ji) - 1, jpi ) + 1 623 ij = ( npb(ji) - 1 ) / jpi + 1 583 624 ! 584 625 WRITE(numout,*) ' alerte 1 ' … … 586 627 WRITE(numout,*) ' heat diffusion in the ice ' 587 628 WRITE(numout,*) ' Category : ', jl 588 WRITE(numout,*) ' zji , zjj : ', zji, zjj589 WRITE(numout,*) ' lat, lon : ', gphit( zji,zjj), glamt(zji,zjj)629 WRITE(numout,*) ' ii , ij : ', ii, ij 630 WRITE(numout,*) ' lat, lon : ', gphit(ii,ij), glamt(ii,ij) 590 631 WRITE(numout,*) ' cons_error : ', cons_error(ji,jl) 591 632 WRITE(numout,*) ' surf_error : ', surf_error(ji,jl) … … 615 656 WRITE(numout,*) ' fc_bo : ', - fc_bo_i (ji) 616 657 WRITE(numout,*) ' foc : ', fbif_1d(ji) 617 WRITE(numout,*) ' fstroc : ', fstroc ( zji,zjj,jl)658 WRITE(numout,*) ' fstroc : ', fstroc (ii,ij,jl) 618 659 WRITE(numout,*) ' i0 : ', i0(ji) 619 660 WRITE(numout,*) ' qsr_ice : ', (1.0-i0(ji))*qsr_ice_1d(ji) … … 651 692 ! 652 693 INTEGER :: ji ! loop indices 653 INTEGER :: zji, zjj, numce ! local integers694 INTEGER :: ii, ij, numce ! local integers 654 695 REAL(wp) :: meance, max_cons_err !local scalar 655 696 !!--------------------------------------------------------------------- … … 669 710 !---------------------------------------- 670 711 DO ji = kideb, kiut 671 zji = MOD( npb(ji) - 1 , jpi ) + 1672 zjj = ( npb(ji) - 1 ) / jpi + 1712 ii = MOD( npb(ji) - 1 , jpi ) + 1 713 ij = ( npb(ji) - 1 ) / jpi + 1 673 714 674 715 fatm (ji,jl) = qnsr_ice_1d(ji) + qsr_ice_1d(ji) ! total heat flux 675 sum_fluxq (ji,jl) = fatm(ji,jl) + fbif_1d(ji) - ftotal_fin(ji) - fstroc( zji,zjj,jl)716 sum_fluxq (ji,jl) = fatm(ji,jl) + fbif_1d(ji) - ftotal_fin(ji) - fstroc(ii,ij,jl) 676 717 cons_error(ji,jl) = ABS( dq_i(ji,jl) * r1_rdtice + sum_fluxq(ji,jl) ) 677 718 END DO … … 706 747 DO ji = kideb, kiut 707 748 IF ( cons_error(ji,jl) .GT. max_cons_err ) THEN 708 zji = MOD( npb(ji) - 1, jpi ) + 1709 zjj = ( npb(ji) - 1 ) / jpi + 1749 ii = MOD( npb(ji) - 1, jpi ) + 1 750 ij = ( npb(ji) - 1 ) / jpi + 1 710 751 ! 711 752 WRITE(numout,*) ' alerte 1 - category : ', jl 712 753 WRITE(numout,*) ' Untolerated conservation error after limthd_ent ' 713 WRITE(numout,*) ' zji , zjj : ', zji, zjj714 WRITE(numout,*) ' lat, lon : ', gphit( zji,zjj), glamt(zji,zjj)754 WRITE(numout,*) ' ii , ij : ', ii, ij 755 WRITE(numout,*) ' lat, lon : ', gphit(ii,ij), glamt(ii,ij) 715 756 WRITE(numout,*) ' * ' 716 757 WRITE(numout,*) ' Ftotal : ', sum_fluxq(ji,jl) … … 724 765 WRITE(numout,*) ' foce : ', fbif_1d(ji) 725 766 WRITE(numout,*) ' fres : ', ftotal_fin(ji) 726 WRITE(numout,*) ' fhbri : ', fhbricat( zji,zjj,jl)767 WRITE(numout,*) ' fhbri : ', fhbricat(ii,ij,jl) 727 768 WRITE(numout,*) ' * ' 728 769 WRITE(numout,*) ' Heat contents --- : ' … … 793 834 INTEGER :: ios ! Local integer output status for namelist read 794 835 NAMELIST/namicethd/ hmelt , hiccrit, fraz_swi, maxfrazb, vfrazb, Cfrazb, & 795 & hicmin, hiclim, amax ,&836 & hicmin, hiclim, & 796 837 & sbeta , parlat, hakspl, hibspl, exld, & 797 838 & hakdif, hnzst , thth , parsub, alphs, betas, & … … 825 866 WRITE(numout,*)' ice thick. corr. to max. energy stored in brine pocket hicmin = ', hicmin 826 867 WRITE(numout,*)' minimum ice thickness hiclim = ', hiclim 827 WRITE(numout,*)' maximum lead fraction amax = ', amax828 868 WRITE(numout,*)' numerical carac. of the scheme for diffusion in ice ' 829 869 WRITE(numout,*)' Cranck-Nicholson (=0.5), implicit (=1), explicit (=0) sbeta = ', sbeta
Note: See TracChangeset
for help on using the changeset viewer.
