- Timestamp:
- 2014-05-27T11:28:12+02:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r4634 r4649 44 44 USE timing ! Timing 45 45 USE cpl_oasis3, ONLY : lk_cpl 46 USE limcons ! conservation tests 46 47 47 48 IMPLICIT NONE … … 86 87 INTEGER :: nbpb ! nb of icy pts for thermo. cal. 87 88 INTEGER :: ii, ij ! temporary dummy loop index 88 REAL(wp) :: zfric_umin = 5e-03_wp ! lower bound for the friction velocity 89 REAL(wp) :: zfric_umax = 2e-02_wp ! upper bound for the friction velocity 90 REAL(wp) :: zinda, zindb, zfric_u ! local scalar 91 REAL(wp) :: zareamin ! - - 92 REAL(wp) :: zchk_v_i, zchk_smv, zchk_e_i, zchk_fs, zchk_fw, zchk_ft, zchk_v_i_b, zchk_smv_b, zchk_e_i_b, zchk_fs_b, zchk_fw_b, zchk_ft_b 93 REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 94 REAL(wp) :: zqld, zqfr 89 REAL(wp) :: zfric_umin = 0._wp ! lower bound for the friction velocity (cice value=5.e-04) 90 REAL(wp) :: zch = 0.0057_wp ! heat transfer coefficient 91 REAL(wp) :: zinda, zindb, zareamin 92 REAL(wp) :: zfric_u, zqld, zqfr 95 93 REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini, zhfx, zqfx 96 94 REAL(wp) :: zhfx_err, ztest 95 ! 96 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 97 97 !!------------------------------------------------------------------- 98 98 IF( nn_timing == 1 ) CALL timing_start('limthd') … … 103 103 zdq(:) = 0._wp ; zq_ini(:) = 0._wp ; zhfx(:) = 0._wp ; zqfx(:) = 0._wp 104 104 105 ! ------------------------------- 106 !- check conservation (C Rousset) 107 IF (ln_limdiahsb) THEN 108 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) 109 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 110 zchk_e_i_b = glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) 111 zchk_fw_b = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) 112 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) 113 zchk_ft_b = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) 114 ENDIF 115 !- check conservation (C Rousset) 116 ! ------------------------------- 105 ! conservation test 106 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 117 107 118 108 !------------------------------------------------------------------------------! … … 158 148 !CDIR NOVERRCHK 159 149 DO ji = 1, jpi 160 zinda = tms(ji,jj) * ( 1. 0- MAX( 0._wp , SIGN( 1._wp , - at_i(ji,jj) + epsi10 ) ) ) ! 0 if no ice150 zinda = tms(ji,jj) * ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - at_i(ji,jj) + epsi10 ) ) ) ! 0 if no ice 161 151 ! 162 152 ! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget … … 165 155 ! ! net downward heat flux from the ice to the ocean, expressed as a function of ocean 166 156 ! ! temperature and turbulent mixing (McPhee, 1992) 167 168 ! friction velocity 169 zfric_u = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin ) 170 171 !-- Energy from the turbulent oceanic heat flux. here the drag will depend on ice thickness and type (0.006) 172 fhtur(ji,jj) = zinda * rau0 * rcp * 0.006 * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2 173 ! clem: why not the following? 174 !fhtur(ji,jj) = zinda * rau0 * rcp * 0.006 * SQRT( ust2s(ji,jj) ) * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) 175 176 !-- Energy received in the lead, zqld is defined everywhere (J.m-2) 177 ! It includes turbulent ocean heat flux (only in the leads, the rest is used for bottom melting) 178 zqld = tms(ji,jj) * rdt_ice * & 179 & ( pfrld(ji,jj) * ( qsr(ji,jj) * oatte(ji,jj) & ! solar heat + clem modif 180 & + qns(ji,jj) & ! non solar heat 181 & + fhtur(ji,jj) ) & ! turbulent ice-ocean heat (0 if no ice) 157 ! 158 ! --- Energy received in the lead, zqld is defined everywhere (J.m-2) --- ! 159 zqld = tms(ji,jj) * rdt_ice * & 160 & ( pfrld(ji,jj) * ( qsr(ji,jj) * oatte(ji,jj) & ! solar heat + clem modif 161 & + qns(ji,jj) ) & ! non solar heat 182 162 ! latent heat of precip (note that precip is included in qns but not in qns_ice) 183 163 & + ( pfrld(ji,jj)**betas - pfrld(ji,jj) ) * sprecip(ji,jj) * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus ) & 184 164 & + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rtt ) ) 185 165 186 !-- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) 166 !-- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- ! 187 167 zqfr = tms(ji,jj) * rau0 * rcp * fse3t_m(ji,jj,1) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) 188 168 … … 192 172 ! If there is ice and leads are warming, then transfer energy from the lead budget and use it for bottom melting 193 173 IF( at_i(ji,jj) > epsi10 .AND. zqld > 0._wp ) THEN 194 fhld (ji,jj) = zqld * r1_rdtice / at_i(ji,jj) ! divided by a _i since this is (re)multiplied by a_i in limthd_dh.F90174 fhld (ji,jj) = zqld * r1_rdtice / at_i(ji,jj) ! divided by at_i since this is (re)multiplied by a_i in limthd_dh.F90 195 175 qlead(ji,jj) = 0._wp 196 176 ENDIF 197 177 ! 198 IF( qlead(ji,jj) == 0._wp ) zqld = 0._wp ; zqfr = 0._wp 199 ! 178 !-- Energy from the turbulent oceanic heat flux --- ! 179 !clem zfric_u = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin ) 180 zfric_u = MAX( SQRT( ust2s(ji,jj) ), zfric_umin ) 181 fhtur(ji,jj) = MAX( 0._wp, zinda * rau0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ) ! W.m-2 182 ! upper bound for fhtur: we do not want SST to drop below Tfreeze. 183 ! So we say that the heat retrieved from the ocean (fhtur+fhld) must be < to the heat necessary to reach Tfreeze (zqfr) 184 ! This is not a clean budget, so that should be corrected at some point 185 fhtur(ji,jj) = zinda * MIN( fhtur(ji,jj), - fhld(ji,jj) - zqfr * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ) 186 200 187 ! ----------------------------------------- 201 188 ! Net heat flux on top of ice-ocean [W.m-2] … … 317 304 CALL tab_2d_1d( nbpb, wfx_sum_1d (1:nbpb), wfx_sum , jpi, jpj, npb(1:nbpb) ) 318 305 CALL tab_2d_1d( nbpb, wfx_sni_1d (1:nbpb), wfx_sni , jpi, jpj, npb(1:nbpb) ) 306 CALL tab_2d_1d( nbpb, wfx_res_1d (1:nbpb), wfx_res , jpi, jpj, npb(1:nbpb) ) 307 CALL tab_2d_1d( nbpb, wfx_spr_1d (1:nbpb), wfx_spr , jpi, jpj, npb(1:nbpb) ) 319 308 320 309 CALL tab_2d_1d( nbpb, sfx_bog_1d (1:nbpb), sfx_bog , jpi, jpj, npb(1:nbpb) ) … … 323 312 CALL tab_2d_1d( nbpb, sfx_sni_1d (1:nbpb), sfx_sni , jpi, jpj, npb(1:nbpb) ) 324 313 CALL tab_2d_1d( nbpb, sfx_bri_1d (1:nbpb), sfx_bri , jpi, jpj, npb(1:nbpb) ) 314 CALL tab_2d_1d( nbpb, sfx_res_1d (1:nbpb), sfx_res , jpi, jpj, npb(1:nbpb) ) 325 315 326 316 CALL tab_2d_1d( nbpb, iatte_1d (1:nbpb), iatte , jpi, jpj, npb(1:nbpb) ) … … 329 319 CALL tab_2d_1d( nbpb, hfx_thd_1d (1:nbpb), hfx_thd , jpi, jpj, npb(1:nbpb) ) 330 320 CALL tab_2d_1d( nbpb, hfx_spr_1d (1:nbpb), hfx_spr , jpi, jpj, npb(1:nbpb) ) 331 CALL tab_2d_1d( nbpb, hfx_tot_1d (1:nbpb), hfx_tot , jpi, jpj, npb(1:nbpb) ) 321 CALL tab_2d_1d( nbpb, hfx_sum_1d (1:nbpb), hfx_sum , jpi, jpj, npb(1:nbpb) ) 322 CALL tab_2d_1d( nbpb, hfx_bom_1d (1:nbpb), hfx_bom , jpi, jpj, npb(1:nbpb) ) 323 CALL tab_2d_1d( nbpb, hfx_bog_1d (1:nbpb), hfx_bog , jpi, jpj, npb(1:nbpb) ) 324 CALL tab_2d_1d( nbpb, hfx_dif_1d (1:nbpb), hfx_dif , jpi, jpj, npb(1:nbpb) ) 325 CALL tab_2d_1d( nbpb, hfx_opw_1d (1:nbpb), hfx_opw , jpi, jpj, npb(1:nbpb) ) 332 326 CALL tab_2d_1d( nbpb, hfx_snw_1d (1:nbpb), hfx_snw , jpi, jpj, npb(1:nbpb) ) 333 327 CALL tab_2d_1d( nbpb, hfx_sub_1d (1:nbpb), hfx_sub , jpi, jpj, npb(1:nbpb) ) … … 365 359 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 366 360 hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_b(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) ) 367 368 361 END DO 369 362 … … 379 372 CALL lim_thd_dh( 1, nbpb, jl ) 380 373 381 ! --- Ice /Snowenthalpy remapping --- !382 CALL lim_thd_ent( 1, nbpb, jl )374 ! --- Ice enthalpy remapping --- ! 375 CALL lim_thd_ent( 1, nbpb, jl, q_i_b(1:nbpb,:) ) 383 376 ! 384 377 ! --- diag error on heat remapping - PART 2 --- ! … … 391 384 392 385 !---------------------------------! 393 ! Ice salinity!386 ! --- Ice salinity --- ! 394 387 !---------------------------------! 395 388 CALL lim_thd_sal( 1, nbpb ) 396 389 397 ! CALL lim_thd_enmelt(1,nbpb) ! computes sea ice energy of melting 390 !---------------------------------! 391 ! --- temperature update --- ! 392 !---------------------------------! 393 CALL lim_thd_temp( 1, nbpb ) 394 398 395 !-------------------------------- 399 396 ! 4.4) Move 1D to 2D vectors … … 424 421 CALL tab_1d_2d( nbpb, wfx_sum , npb, wfx_sum_1d(1:nbpb) , jpi, jpj ) 425 422 CALL tab_1d_2d( nbpb, wfx_sni , npb, wfx_sni_1d(1:nbpb) , jpi, jpj ) 423 CALL tab_1d_2d( nbpb, wfx_res , npb, wfx_res_1d(1:nbpb) , jpi, jpj ) 424 CALL tab_1d_2d( nbpb, wfx_spr , npb, wfx_spr_1d(1:nbpb) , jpi, jpj ) 426 425 427 426 CALL tab_1d_2d( nbpb, sfx_bog , npb, sfx_bog_1d(1:nbpb) , jpi, jpj ) … … 429 428 CALL tab_1d_2d( nbpb, sfx_sum , npb, sfx_sum_1d(1:nbpb) , jpi, jpj ) 430 429 CALL tab_1d_2d( nbpb, sfx_sni , npb, sfx_sni_1d(1:nbpb) , jpi, jpj ) 430 CALL tab_1d_2d( nbpb, sfx_res , npb, sfx_res_1d(1:nbpb) , jpi, jpj ) 431 431 ! 432 432 IF( num_sal == 2 ) THEN … … 436 436 CALL tab_1d_2d( nbpb, hfx_thd , npb, hfx_thd_1d(1:nbpb) , jpi, jpj ) 437 437 CALL tab_1d_2d( nbpb, hfx_spr , npb, hfx_spr_1d(1:nbpb) , jpi, jpj ) 438 CALL tab_1d_2d( nbpb, hfx_tot , npb, hfx_tot_1d(1:nbpb) , jpi, jpj ) 438 CALL tab_1d_2d( nbpb, hfx_sum , npb, hfx_sum_1d(1:nbpb) , jpi, jpj ) 439 CALL tab_1d_2d( nbpb, hfx_bom , npb, hfx_bom_1d(1:nbpb) , jpi, jpj ) 440 CALL tab_1d_2d( nbpb, hfx_bog , npb, hfx_bog_1d(1:nbpb) , jpi, jpj ) 441 CALL tab_1d_2d( nbpb, hfx_dif , npb, hfx_dif_1d(1:nbpb) , jpi, jpj ) 442 CALL tab_1d_2d( nbpb, hfx_opw , npb, hfx_opw_1d(1:nbpb) , jpi, jpj ) 439 443 CALL tab_1d_2d( nbpb, hfx_snw , npb, hfx_snw_1d(1:nbpb) , jpi, jpj ) 440 444 CALL tab_1d_2d( nbpb, hfx_sub , npb, hfx_sub_1d(1:nbpb) , jpi, jpj ) … … 521 525 ENDIF 522 526 ! 523 ! ------------------------------- 524 !- check conservation (C Rousset) 525 IF (ln_limdiahsb) THEN 526 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 527 zchk_fw = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fw_b 528 zchk_ft = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) - zchk_ft_b 529 530 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b ) * r1_rdtice - zchk_fw 531 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 532 zchk_e_i = glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) * r1_rdtice - zchk_e_i_b * r1_rdtice + zchk_ft 533 534 zchk_vmin = glob_min(v_i) 535 zchk_amax = glob_max(SUM(a_i,dim=3)) 536 zchk_amin = glob_min(a_i) 537 538 IF(lwp) THEN 539 IF ( ABS( zchk_v_i ) > 1.e-4 ) WRITE(numout,*) 'violation volume [kg/day] (limthd) = ',(zchk_v_i * rday) 540 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limthd) = ',(zchk_smv * rday) 541 IF ( ABS( zchk_e_i ) > 1.e-2 ) WRITE(numout,*) 'violation enthalpy [1e9 J] (limthd) = ',(zchk_e_i) 542 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limthd) = ',(zchk_vmin * 1.e-3) 543 IF ( zchk_amax > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limthd) = ',zchk_amax 544 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limthd) = ',zchk_amin 545 ENDIF 546 ENDIF 547 !- check conservation (C Rousset) 548 ! ------------------------------- 527 ! conservation test 528 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 549 529 ! 550 530 CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx, zqfx ) … … 558 538 !! *** ROUTINE lim_thd_enmelt *** 559 539 !! 560 !! ** Purpose : Computes sea ice energy of melting q_i (J.m-3) 540 !! ** Purpose : Computes sea ice energy of melting q_i (J.m-3) from temperature 561 541 !! 562 542 !! ** Method : Formula (Bitz and Lipscomb, 1999) … … 565 545 !! 566 546 INTEGER :: ji, jk ! dummy loop indices 567 REAL(wp) :: ztmelts ! local scalar547 REAL(wp) :: ztmelts, zindb ! local scalar 568 548 !!------------------------------------------------------------------- 569 549 ! 570 550 DO jk = 1, nlay_i ! Sea ice energy of melting 571 551 DO ji = kideb, kiut 572 ztmelts = - tmut * s_i_b(ji,jk) + rtt 573 q_i_b(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_b(ji,jk) ) & 574 & + lfus * ( 1.0 - (ztmelts-rtt) / MIN( t_i_b(ji,jk)-rtt, -epsi10 ) ) & 575 & - rcp * ( ztmelts-rtt ) ) 552 ztmelts = - tmut * s_i_b(ji,jk) + rtt 553 zindb = MAX( 0._wp , SIGN( 1._wp , -(t_i_b(ji,jk) - rtt) - epsi10 ) ) 554 q_i_b(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_b(ji,jk) ) & 555 & + lfus * ( 1.0 - zindb * ( ztmelts-rtt ) / MIN( t_i_b(ji,jk)-rtt, -epsi10 ) ) & 556 & - rcp * ( ztmelts-rtt ) ) 576 557 END DO 577 558 END DO … … 584 565 END SUBROUTINE lim_thd_enmelt 585 566 567 SUBROUTINE lim_thd_temp( kideb, kiut ) 568 !!----------------------------------------------------------------------- 569 !! *** ROUTINE lim_thd_temp *** 570 !! 571 !! ** Purpose : Computes sea ice temperature (Kelvin) from enthalpy 572 !! 573 !! ** Method : Formula (Bitz and Lipscomb, 1999) 574 !!------------------------------------------------------------------- 575 INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop 576 !! 577 INTEGER :: ji, jk ! dummy loop indices 578 REAL(wp) :: ztmelts, zswitch, zaaa, zbbb, zccc, zdiscrim ! local scalar 579 !!------------------------------------------------------------------- 580 ! Recover ice temperature 581 DO jk = 1, nlay_i 582 DO ji = kideb, kiut 583 ztmelts = -tmut * s_i_b(ji,jk) + rtt 584 ! Conversion q(S,T) -> T (second order equation) 585 zaaa = cpic 586 zbbb = ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_b(ji,jk) / rhoic - lfus 587 zccc = lfus * ( ztmelts - rtt ) 588 zdiscrim = SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) ) 589 t_i_b(ji,jk) = rtt - ( zbbb + zdiscrim ) / ( 2._wp * zaaa ) 590 591 ! mask temperature 592 zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_b(ji) ) ) 593 t_i_b(ji,jk) = zswitch * t_i_b(ji,jk) + ( 1._wp - zswitch ) * rtt 594 END DO 595 END DO 596 597 END SUBROUTINE lim_thd_temp 586 598 587 599 SUBROUTINE lim_thd_init
Note: See TracChangeset
for help on using the changeset viewer.
