- Timestamp:
- 2017-10-04T09:19:23+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/RK3_SRC/SBC/sbcblk.F90
r7753 r8586 40 40 USE lib_fortran ! to use key_nosignedzero 41 41 #if defined key_lim3 42 USE ice , ONLY : u_ice, v_ice, jpl, pfrld, a_i_b, at_i_b 43 USE limthd_dh ! for CALL lim_thd_snwblow 44 #elif defined key_lim2 45 USE ice_2 , ONLY : u_ice, v_ice 46 USE par_ice_2 ! LIM-2 parameters 42 USE ice , ONLY : u_ice, v_ice, jpl, a_i_b, at_i_b 43 USE icethd_dh ! for CALL ice_thd_snwblow 47 44 #endif 48 45 USE sbcblk_algo_ncar ! => turb_ncar : NCAR - CORE (Large & Yeager, 2009) … … 54 51 USE in_out_manager ! I/O manager 55 52 USE lib_mpp ! distribued memory computing library 56 USE wrk_nemo ! work arrays57 53 USE timing ! Timing 58 54 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 64 60 PUBLIC sbc_blk_init ! called in sbcmod 65 61 PUBLIC sbc_blk ! called in sbcmod 66 #if defined key_lim 2 || defined key_lim367 PUBLIC blk_ice_tau ! routine called in sbc_ice_limmodule68 PUBLIC blk_ice_flx ! routine called in sbc_ice_limmodule62 #if defined key_lim3 63 PUBLIC blk_ice_tau ! routine called in icestp module 64 PUBLIC blk_ice_flx ! routine called in icestp module 69 65 #endif 70 66 … … 111 107 REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements 112 108 REAL(wp) :: rn_zu ! z(u) : height of wind measurements 109 !!gm ref namelist initialize it so remove the setting to false below 113 110 LOGICAL :: ln_Cd_L12 = .FALSE. ! Modify the drag ice-atm and oce-atm depending on ice concentration (from Lupkes et al. JGR2012) 114 111 ! … … 360 357 INTEGER :: ji, jj ! dummy loop indices 361 358 REAL(wp) :: zztmp ! local variable 362 REAL(wp), DIMENSION(:,:), POINTER :: zwnd_i, zwnd_j ! wind speed components at T-point 363 REAL(wp), DIMENSION(:,:), POINTER :: zsq ! specific humidity at pst 364 REAL(wp), DIMENSION(:,:), POINTER :: zqlw, zqsb ! long wave and sensible heat fluxes 365 REAL(wp), DIMENSION(:,:), POINTER :: zqla, zevap ! latent heat fluxes and evaporation 366 REAL(wp), DIMENSION(:,:), POINTER :: Cd ! transfer coefficient for momentum (tau) 367 REAL(wp), DIMENSION(:,:), POINTER :: Ch ! transfer coefficient for sensible heat (Q_sens) 368 REAL(wp), DIMENSION(:,:), POINTER :: Ce ! tansfert coefficient for evaporation (Q_lat) 369 REAL(wp), DIMENSION(:,:), POINTER :: zst ! surface temperature in Kelvin 370 REAL(wp), DIMENSION(:,:), POINTER :: zt_zu ! air temperature at wind speed height 371 REAL(wp), DIMENSION(:,:), POINTER :: zq_zu ! air spec. hum. at wind speed height 372 REAL(wp), DIMENSION(:,:), POINTER :: zU_zu ! bulk wind speed at height zu [m/s] 373 REAL(wp), DIMENSION(:,:), POINTER :: ztpot ! potential temperature of air at z=rn_zqt [K] 374 REAL(wp), DIMENSION(:,:), POINTER :: zrhoa ! density of air [kg/m^3] 375 !!--------------------------------------------------------------------- 376 ! 377 IF( nn_timing == 1 ) CALL timing_start('blk_oce') 378 ! 379 CALL wrk_alloc( jpi,jpj, zwnd_i, zwnd_j, zsq, zqlw, zqsb, zqla, zevap ) 380 CALL wrk_alloc( jpi,jpj, Cd, Ch, Ce, zst, zt_zu, zq_zu ) 381 CALL wrk_alloc( jpi,jpj, zU_zu, ztpot, zrhoa ) 382 ! 383 359 REAL(wp), DIMENSION(jpi,jpj) :: zwnd_i, zwnd_j ! wind speed components at T-point 360 REAL(wp), DIMENSION(jpi,jpj) :: zsq ! specific humidity at pst 361 REAL(wp), DIMENSION(jpi,jpj) :: zqlw, zqsb ! long wave and sensible heat fluxes 362 REAL(wp), DIMENSION(jpi,jpj) :: zqla, zevap ! latent heat fluxes and evaporation 363 REAL(wp), DIMENSION(jpi,jpj) :: zCd ! transfer coefficient for momentum (tau) 364 REAL(wp), DIMENSION(jpi,jpj) :: zCh ! transfer coefficient for sensible heat (Q_sens) 365 REAL(wp), DIMENSION(jpi,jpj) :: zCe ! tansfert coefficient for evaporation (Q_lat) 366 REAL(wp), DIMENSION(jpi,jpj) :: zst ! surface temperature in Kelvin 367 REAL(wp), DIMENSION(jpi,jpj) :: zt_zu ! air temperature at wind speed height 368 REAL(wp), DIMENSION(jpi,jpj) :: zq_zu ! air spec. hum. at wind speed height 369 REAL(wp), DIMENSION(jpi,jpj) :: zU_zu ! bulk wind speed at height zu [m/s] 370 REAL(wp), DIMENSION(jpi,jpj) :: ztpot ! potential temperature of air at z=rn_zqt [K] 371 REAL(wp), DIMENSION(jpi,jpj) :: zrhoa ! density of air [kg/m^3] 372 !!--------------------------------------------------------------------- 373 ! 374 IF( ln_timing ) CALL timing_start('blk_oce') 375 ! 384 376 ! local scalars ( place there for vector optimisation purposes) 385 377 zst(:,:) = pst(:,:) + rt0 ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) … … 443 435 ! 444 436 CASE( np_NCAR ) ; CALL turb_ncar ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! NCAR-COREv2 445 & Cd, Ch,Ce, zt_zu, zq_zu, zU_zu )437 & zCd, zCh, zCe, zt_zu, zq_zu, zU_zu ) 446 438 CASE( np_COARE_3p0 ) ; CALL turb_coare ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! COARE v3.0 447 & Cd, Ch,Ce, zt_zu, zq_zu, zU_zu )439 & zCd, zCh, zCe, zt_zu, zq_zu, zU_zu ) 448 440 CASE( np_COARE_3p5 ) ; CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! COARE v3.5 449 & Cd, Ch,Ce, zt_zu, zq_zu, zU_zu )441 & zCd, zCh, zCe, zt_zu, zq_zu, zU_zu ) 450 442 CASE( np_ECMWF ) ; CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! ECMWF 451 & Cd, Ch,Ce, zt_zu, zq_zu, zU_zu )443 & zCd, zCh, zCe, zt_zu, zq_zu, zU_zu ) 452 444 CASE DEFAULT 453 445 CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) … … 461 453 END IF 462 454 463 Cd_oce(:,:) = Cd(:,:)! record value of pure ocean-atm. drag (clem)455 Cd_oce(:,:) = zCd(:,:) ! record value of pure ocean-atm. drag (clem) 464 456 465 457 DO jj = 1, jpj ! tau module, i and j component 466 458 DO ji = 1, jpi 467 zztmp = zrhoa(ji,jj) * zU_zu(ji,jj) * Cd(ji,jj) ! using bulk wind speed459 zztmp = zrhoa(ji,jj) * zU_zu(ji,jj) * zCd(ji,jj) ! using bulk wind speed 468 460 taum (ji,jj) = zztmp * wndm (ji,jj) 469 461 zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) … … 500 492 IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN 501 493 !! q_air and t_air are given at 10m (wind reference height) 502 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)* Ce(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) )! Evaporation, using bulk wind speed503 zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)* Ch(:,:)*(zst(:,:) - ztpot(:,:)) ! Sensible Heat, using bulk wind speed494 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)*zCe(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed 495 zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*zCh(:,:)*(zst(:,:) - ztpot(:,:) ) ! Sensible Heat, using bulk wind speed 504 496 ELSE 505 497 !! q_air and t_air are not given at 10m (wind reference height) 506 498 ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!! 507 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)* Ce(:,:)*(zsq(:,:) - zq_zu(:,:) ) ) ! Evaporation ! using bulk wind speed508 zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)* Ch(:,:)*(zst(:,:) - zt_zu(:,:) ) ! Sensible Heat ! using bulk wind speed499 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)*zCe(:,:)*(zsq(:,:) - zq_zu(:,:) ) ) ! Evaporation ! using bulk wind speed 500 zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*zCh(:,:)*(zst(:,:) - zt_zu(:,:) ) ! Sensible Heat ! using bulk wind speed 509 501 ENDIF 510 502 … … 513 505 514 506 IF(ln_ctl) THEN 515 CALL prt_ctl( tab2d_1=zqla , clinfo1=' blk_oce: zqla : ', tab2d_2= Ce , clinfo2=' Ce : ' )516 CALL prt_ctl( tab2d_1=zqsb , clinfo1=' blk_oce: zqsb : ', tab2d_2= Ch , clinfo2=' Ch : ' )507 CALL prt_ctl( tab2d_1=zqla , clinfo1=' blk_oce: zqla : ', tab2d_2=zCe , clinfo2=' Ce : ' ) 508 CALL prt_ctl( tab2d_1=zqsb , clinfo1=' blk_oce: zqsb : ', tab2d_2=zCh , clinfo2=' Ch : ' ) 517 509 CALL prt_ctl( tab2d_1=zqlw , clinfo1=' blk_oce: zqlw : ', tab2d_2=qsr, clinfo2=' qsr : ' ) 518 510 CALL prt_ctl( tab2d_1=zsq , clinfo1=' blk_oce: zsq : ', tab2d_2=zst, clinfo2=' zst : ' ) … … 565 557 ENDIF 566 558 ! 567 CALL wrk_dealloc( jpi,jpj, zwnd_i, zwnd_j, zsq, zqlw, zqsb, zqla, zevap ) 568 CALL wrk_dealloc( jpi,jpj, Cd, Ch, Ce, zst, zt_zu, zq_zu ) 569 CALL wrk_dealloc( jpi,jpj, zU_zu, ztpot, zrhoa ) 570 ! 571 IF( nn_timing == 1 ) CALL timing_stop('blk_oce') 559 IF( ln_timing ) CALL timing_stop('blk_oce') 572 560 ! 573 561 END SUBROUTINE blk_oce 574 562 575 #if defined key_lim 2 || defined key_lim3563 #if defined key_lim3 576 564 577 565 SUBROUTINE blk_ice_tau … … 586 574 !!--------------------------------------------------------------------- 587 575 INTEGER :: ji, jj ! dummy loop indices 588 ! 589 REAL(wp), DIMENSION(:,:) , POINTER :: zrhoa 590 ! 591 REAL(wp) :: zwnorm_f, zwndi_f , zwndj_f ! relative wind module and components at F-point 592 REAL(wp) :: zwndi_t , zwndj_t ! relative wind components at T-point 593 REAL(wp), DIMENSION(:,:), POINTER :: Cd ! transfer coefficient for momentum (tau) 594 !!--------------------------------------------------------------------- 595 ! 596 IF( nn_timing == 1 ) CALL timing_start('blk_ice_tau') 597 ! 598 CALL wrk_alloc( jpi,jpj, zrhoa ) 599 CALL wrk_alloc( jpi,jpj, Cd ) 600 601 Cd(:,:) = Cd_ice 602 603 ! Make ice-atm. drag dependent on ice concentration (see Lupkes et al. 2012) (clem) 604 #if defined key_lim3 576 REAL(wp) :: zwndi_f , zwndj_f, zwnorm_f ! relative wind module and components at F-point 577 REAL(wp) :: zwndi_t , zwndj_t ! relative wind components at T-point 578 REAL(wp), DIMENSION(jpi,jpj) :: zCd, zrhoa ! transfer coefficient for momentum (tau) 579 !!--------------------------------------------------------------------- 580 ! 581 IF( ln_timing ) CALL timing_start('blk_ice_tau') 582 ! 605 583 IF( ln_Cd_L12 ) THEN 606 CALL Cdn10_Lupkes2012( Cd ) ! calculate new drag from Lupkes(2012) equations 607 ENDIF 608 #endif 584 CALL Cdn10_Lupkes2012( zCd ) ! air-ice drag = F(ice concentration) (see Lupkes et al., 2012) 585 ELSE 586 zCd(:,:) = Cd_ice ! constant air-ice drag 587 ENDIF 609 588 610 589 ! local scalars ( place there for vector optimisation purposes) … … 632 611 zwndj_f = 0.25 * ( sf(jp_wndj)%fnow(ji-1,jj ,1) + sf(jp_wndj)%fnow(ji ,jj ,1) & 633 612 & + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji ,jj-1,1) ) - rn_vfac * v_ice(ji,jj) 634 zwnorm_f = zrhoa(ji,jj) * Cd(ji,jj) * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f )613 zwnorm_f = zrhoa(ji,jj) * zCd(ji,jj) * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 635 614 ! ... ice stress at I-point 636 615 utau_ice(ji,jj) = zwnorm_f * zwndi_f … … 658 637 DO jj = 2, jpjm1 659 638 DO ji = fs_2, fs_jpim1 ! vect. opt. 660 utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd(ji,jj) * ( wndm_ice(ji+1,jj ) + wndm_ice(ji,jj) ) &639 utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * zCd(ji,jj) * ( wndm_ice(ji+1,jj ) + wndm_ice(ji,jj) ) & 661 640 & * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) ) 662 vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd(ji,jj) * ( wndm_ice(ji,jj+1 ) + wndm_ice(ji,jj) ) &641 vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * zCd(ji,jj) * ( wndm_ice(ji,jj+1 ) + wndm_ice(ji,jj) ) & 663 642 & * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) ) 664 643 END DO … … 669 648 ! 670 649 END SELECT 671 650 ! 672 651 IF(ln_ctl) THEN 673 652 CALL prt_ctl(tab2d_1=utau_ice , clinfo1=' blk_ice: utau_ice : ', tab2d_2=vtau_ice , clinfo2=' vtau_ice : ') 674 653 CALL prt_ctl(tab2d_1=wndm_ice , clinfo1=' blk_ice: wndm_ice : ') 675 654 ENDIF 676 677 IF( nn_timing == 1 )CALL timing_stop('blk_ice_tau')678 655 ! 656 IF( ln_timing ) CALL timing_stop('blk_ice_tau') 657 ! 679 658 END SUBROUTINE blk_ice_tau 680 659 … … 699 678 REAL(wp) :: zcoef_dqlw, zcoef_dqla ! - - 700 679 REAL(wp) :: zztmp, z1_lsub ! - - 701 REAL(wp), DIMENSION(:,:,:), POINTER :: z_qlw ! long wave heat flux over ice 702 REAL(wp), DIMENSION(:,:,:), POINTER :: z_qsb ! sensible heat flux over ice 703 REAL(wp), DIMENSION(:,:,:), POINTER :: z_dqlw ! long wave heat sensitivity over ice 704 REAL(wp), DIMENSION(:,:,:), POINTER :: z_dqsb ! sensible heat sensitivity over ice 705 REAL(wp), DIMENSION(:,:) , POINTER :: zevap, zsnw ! evaporation and snw distribution after wind blowing (LIM3) 706 REAL(wp), DIMENSION(:,:) , POINTER :: zrhoa 707 REAL(wp), DIMENSION(:,:) , POINTER :: Cd ! transfer coefficient for momentum (tau) 708 !!--------------------------------------------------------------------- 709 ! 710 IF( nn_timing == 1 ) CALL timing_start('blk_ice_flx') 711 ! 712 CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb ) 713 CALL wrk_alloc( jpi,jpj, zrhoa) 714 CALL wrk_alloc( jpi,jpj, Cd ) 715 716 Cd(:,:) = Cd_ice 717 718 ! Make ice-atm. drag dependent on ice concentration (see Lupkes et al. 2012) (clem) 719 #if defined key_lim3 680 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qlw ! long wave heat flux over ice 681 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qsb ! sensible heat flux over ice 682 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqlw ! long wave heat sensitivity over ice 683 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqsb ! sensible heat sensitivity over ice 684 REAL(wp), DIMENSION(jpi,jpj) :: zevap, zsnw ! evaporation and snw distribution after wind blowing (LIM3) 685 REAL(wp), DIMENSION(jpi,jpj) :: zrhoa 686 REAL(wp), DIMENSION(jpi,jpj) :: zCd ! transfer coefficient for momentum (tau) 687 !!--------------------------------------------------------------------- 688 ! 689 IF( ln_timing ) CALL timing_start('blk_ice_flx') 690 ! 720 691 IF( ln_Cd_L12 ) THEN 721 CALL Cdn10_Lupkes2012( Cd ) ! calculate new drag from Lupkes(2012) equations 722 ENDIF 723 #endif 692 CALL Cdn10_Lupkes2012( zCd ) ! air-ice drag = F(ice concentration) (see Lupkes et al., 2012) 693 ELSE 694 zCd(:,:) = Cd_ice ! constant air-ice drag 695 ENDIF 724 696 725 697 ! … … 754 726 ! ... turbulent heat fluxes 755 727 ! Sensible Heat 756 z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Cd(ji,jj) * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) )728 z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * zCd(ji,jj) * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) ) 757 729 ! Latent Heat 758 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls * Cd(ji,jj) * wndm_ice(ji,jj) &730 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls * zCd(ji,jj) * wndm_ice(ji,jj) & 759 731 & * ( 11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) ) 760 732 ! Latent heat sensitivity for ice (Dqla/Dt) 761 733 IF( qla_ice(ji,jj,jl) > 0._wp ) THEN 762 dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Cd(ji,jj) * wndm_ice(ji,jj) / ( zst2 ) * EXP( -5897.8 / ptsu(ji,jj,jl) )734 dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * zCd(ji,jj) * wndm_ice(ji,jj) / ( zst2 ) * EXP( -5897.8 / ptsu(ji,jj,jl) ) 763 735 ELSE 764 736 dqla_ice(ji,jj,jl) = 0._wp … … 766 738 767 739 ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) 768 z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Cd(ji,jj) * wndm_ice(ji,jj)740 z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * zCd(ji,jj) * wndm_ice(ji,jj) 769 741 770 742 ! ----------------------------! … … 786 758 CALL iom_put( 'precip' , tprecip * 86400. ) ! Total precipitation 787 759 788 #if defined key_lim3789 CALL wrk_alloc( jpi,jpj, zevap, zsnw )790 791 760 ! --- evaporation --- ! 792 761 z1_lsub = 1._wp / Lsub … … 797 766 ! --- evaporation minus precipitation --- ! 798 767 zsnw(:,:) = 0._wp 799 CALL lim_thd_snwblow( pfrld, zsnw ) ! snow distribution over ice after wind blowing800 emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )768 CALL ice_thd_snwblow( (1.-at_i_b(:,:)), zsnw ) ! snow distribution over ice after wind blowing 769 emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw ) 801 770 emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw 802 771 emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:) 803 772 804 773 ! --- heat flux associated with emp --- ! 805 qemp_oce(:,:) = - pfrld(:,:) * zevap(:,:) * sst_m(:,:) * rcp& ! evap at sst774 qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp & ! evap at sst 806 775 & + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp & ! liquid precip at Tair 807 776 & + sprecip(:,:) * ( 1._wp - zsnw ) * & ! solid precip at min(Tair,Tsnow) … … 811 780 812 781 ! --- total solar and non solar fluxes --- ! 813 qns_tot(:,:) = pfrld(:,:) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) + qemp_ice(:,:) + qemp_oce(:,:) 814 qsr_tot(:,:) = pfrld(:,:) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 ) 782 qns_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) & 783 & + qemp_ice(:,:) + qemp_oce(:,:) 784 qsr_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 ) 815 785 816 786 ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- ! … … 820 790 DO jl = 1, jpl 821 791 qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) ) 822 792 ! ! But we do not have Tice => consider it at 0degC => evap=0 823 793 END DO 824 825 CALL wrk_dealloc( jpi,jpj, zevap, zsnw )826 #endif827 794 828 795 !-------------------------------------------------------------------- … … 833 800 fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 834 801 fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 835 !836 802 ! 837 803 IF(ln_ctl) THEN … … 843 809 CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip : ', tab2d_2=sprecip , clinfo2=' sprecip : ') 844 810 ENDIF 845 846 CALL wrk_dealloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb ) 847 CALL wrk_dealloc( jpi,jpj, zrhoa ) 848 CALL wrk_dealloc( jpi,jpj, Cd ) 849 ! 850 IF( nn_timing == 1 ) CALL timing_stop('blk_ice_flx') 851 811 ! 812 IF( ln_timing ) CALL timing_stop('blk_ice_flx') 813 ! 852 814 END SUBROUTINE blk_ice_flx 853 815 … … 971 933 END FUNCTION L_vap 972 934 973 974 935 #if defined key_lim3 975 SUBROUTINE Cdn10_Lupkes2012( Cd ) 936 937 SUBROUTINE Cdn10_Lupkes2012( pCd ) 976 938 !!---------------------------------------------------------------------- 977 939 !! *** ROUTINE Cdn10_Lupkes2012 *** … … 1003 965 !! 1004 966 !!---------------------------------------------------------------------- 1005 REAL(wp), DIMENSION(:,:), INTENT( inout) :: Cd967 REAL(wp), DIMENSION(:,:), INTENT( out) :: pCd ! air-ice drag coefficient 1006 968 REAL(wp), PARAMETER :: zCe = 2.23e-03_wp 1007 969 REAL(wp), PARAMETER :: znu = 1._wp … … 1011 973 !!---------------------------------------------------------------------- 1012 974 zcoef = znu + 1._wp / ( 10._wp * zbeta ) 1013 975 ! 1014 976 ! generic drag over a cell partly covered by ice 1015 !! Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) + & ! pure ocean drag1016 !! & Cd_ice * at_i_b(:,:) + & ! pure ice drag1017 !! & zCe * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**zmu ! change due to sea-ice morphology977 !!pCd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) + & ! pure ocean drag 978 !! & Cd_ice * at_i_b(:,:) + & ! pure ice drag 979 !! & zCe * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**zmu ! change due to sea-ice morphology 1018 980 1019 981 ! ice-atm drag 1020 Cd(:,:) = Cd_ice + & ! pure ice drag1021 & zCe * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp) ! change due to sea-ice morphology1022 982 pCd(:,:) = Cd_ice + & ! pure ice drag 983 & zCe * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp) ! change due to sea-ice morphology 984 ! 1023 985 END SUBROUTINE Cdn10_Lupkes2012 1024 986 #endif
Note: See TracChangeset
for help on using the changeset viewer.
