New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 11111 for NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90 – NEMO

Ignore:
Timestamp:
2019-06-14T15:55:28+02:00 (5 years ago)
Author:
laurent
Message:

LB: skin temperature now used by ECMWF bulk algorithm, and can be xios-ed ; new module "SBC/sbcblk_skin.F90" ; all physical constants defined in SBC now moved to "DOM/phycst.F90"; all air-properties and other ABL functions gathered in a new module: "SBC/sbcblk_phymbl.F90".

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90

    r10535 r11111  
    1515   !!            3.7  !  2014-06  (L. Brodeau)  simplification and optimization of CORE bulk 
    1616   !!            4.0  !  2016-06  (L. Brodeau)  sbcblk_core becomes sbcblk and is not restricted to the CORE algorithm anymore 
    17    !!                 !                        ==> based on AeroBulk (http://aerobulk.sourceforge.net/) 
     17   !!                 !                        ==> based on AeroBulk (https://github.com/brodeau/aerobulk/) 
    1818   !!            4.0  !  2016-10  (G. Madec)  introduce a sbc_blk_init routine 
    1919   !!            4.0  !  2016-10  (M. Vancoppenolle)  Introduce conduction flux emulator (M. Vancoppenolle)  
     
    2424   !!   sbc_blk       : bulk formulation as ocean surface boundary condition 
    2525   !!   blk_oce       : computes momentum, heat and freshwater fluxes over ocean 
    26    !!   rho_air       : density of (moist) air (depends on T_air, q_air and SLP 
    27    !!   cp_air        : specific heat of (moist) air (depends spec. hum. q_air) 
    28    !!   q_sat         : saturation humidity as a function of SLP and temperature 
    29    !!   L_vap         : latent heat of vaporization of water as a function of temperature 
    3026   !!             sea-ice case only :  
    3127   !!   blk_ice_tau   : provide the air-ice stress 
     
    6056   USE prtctl         ! Print control 
    6157 
     58   USE sbcblk_phymbl  !LB: all thermodynamics functions, rho_air, q_sat, etc... 
     59 
     60    
    6261   IMPLICIT NONE 
    6362   PRIVATE 
     
    7069   PUBLIC   blk_ice_qcn   ! routine called in icesbc 
    7170#endif  
    72  
    73 !!Lolo: should ultimately be moved in the module with all physical constants ? 
    74 !!gm  : In principle, yes. 
    75    REAL(wp), PARAMETER ::   Cp_dry = 1005.0       !: Specic heat of dry air, constant pressure      [J/K/kg] 
    76    REAL(wp), PARAMETER ::   Cp_vap = 1860.0       !: Specic heat of water vapor, constant pressure  [J/K/kg] 
    77    REAL(wp), PARAMETER ::   R_dry = 287.05_wp     !: Specific gas constant for dry air              [J/K/kg] 
    78    REAL(wp), PARAMETER ::   R_vap = 461.495_wp    !: Specific gas constant for water vapor          [J/K/kg] 
    79    REAL(wp), PARAMETER ::   reps0 = R_dry/R_vap   !: ratio of gas constant for dry air and water vapor => ~ 0.622 
    80    REAL(wp), PARAMETER ::   rctv0 = R_vap/R_dry   !: for virtual temperature (== (1-eps)/eps) => ~ 0.608 
    8171 
    8272   INTEGER , PARAMETER ::   jpfld   =10           ! maximum number of files to read 
     
    9484   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input fields (file informations, fields read) 
    9585 
    96    !                                             !!! Bulk parameters 
    97    REAL(wp), PARAMETER ::   cpa    = 1000.5         ! specific heat of air (only used for ice fluxes now...) 
    98    REAL(wp), PARAMETER ::   Ls     =    2.839e6     ! latent heat of sublimation 
    99    REAL(wp), PARAMETER ::   Stef   =    5.67e-8     ! Stefan Boltzmann constant 
    100    REAL(wp), PARAMETER ::   Cd_ice =    1.4e-3      ! transfer coefficient over ice 
    101    REAL(wp), PARAMETER ::   albo   =    0.066       ! ocean albedo assumed to be constant 
    102    ! 
    10386   !                           !!* Namelist namsbc_blk : bulk parameters 
    10487   LOGICAL  ::   ln_NCAR        ! "NCAR"      algorithm   (Large and Yeager 2008) 
     
    124107   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cdn_oce, chn_oce, cen_oce ! needed by Lupkes 2015 bulk scheme 
    125108 
     109   !LB: 
     110   LOGICAL  ::   ln_skin        ! use the cool-skin / warm-layer parameterization (only available in ECMWF and COARE algorithms) !LB 
     111   !LB: => sbc_oce.F90: REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tsk                       ! sea-surface skin temperature out of the cool-skin/warm-layer parameterization [Celsius] 
     112   !LB. 
     113    
    126114   INTEGER  ::   nblk           ! choice of the bulk algorithm 
    127115   !                            ! associated indices: 
     
    150138      IF( sbc_blk_alloc /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_alloc: failed to allocate arrays' ) 
    151139   END FUNCTION sbc_blk_alloc 
     140 
     141   !LB: 
     142   INTEGER FUNCTION sbc_blk_cswl_alloc() 
     143      !!------------------------------------------------------------------- 
     144      !!             ***  ROUTINE sbc_blk_cswl_alloc *** 
     145      !!------------------------------------------------------------------- 
     146      !PRINT *, '*** LB: allocating tsk!' 
     147      ALLOCATE( tsk(jpi,jpj), STAT=sbc_blk_cswl_alloc ) 
     148      !PRINT *, '*** LB: done!' 
     149      CALL mpp_sum ( 'sbcblk', sbc_blk_cswl_alloc ) 
     150      IF( sbc_blk_cswl_alloc /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_cswl_alloc: failed to allocate arrays' ) 
     151   END FUNCTION sbc_blk_cswl_alloc 
     152   !LB. 
    152153 
    153154 
     
    173174         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF,             &   ! bulk algorithm 
    174175         &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         &  
    175          &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15 
     176         &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15,           & 
     177         &                 ln_skin                                                        ! cool-skin / warm-layer param. !LB 
    176178      !!--------------------------------------------------------------------- 
    177179      ! 
     
    198200      IF( ln_ECMWF     ) THEN   ;   nblk =  np_ECMWF       ;   ioptio = ioptio + 1   ;   ENDIF 
    199201      ! 
     202      !LB: 
     203      !                             !** initialization of the cool-skin / warm-layer parametrization 
     204      IF( ln_skin ) THEN 
     205         IF ( ln_NCAR ) CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with NCAR algorithm!' ) 
     206         !                       ! allocate array(s) for cool-skin/warm-layer param. 
     207         IF( sbc_blk_cswl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' ) 
     208      END IF 
     209      !LB. 
     210          
    200211      IF( ioptio /= 1 )   CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' ) 
    201212      ! 
     
    279290         END SELECT 
    280291         ! 
     292         WRITE(numout,*) 
     293         WRITE(numout,*) '      use cool-skin/warm-layer parameterization (SSST)     ln_skin     = ', ln_skin !LB 
     294         ! 
    281295      ENDIF 
    282296      ! 
     
    413427 
    414428      ! ----------------------------------------------------------------------------- ! 
    415       !      I   Radiative FLUXES                                                     ! 
     429      !      I   Solar FLUX                                                           ! 
    416430      ! ----------------------------------------------------------------------------- ! 
    417431 
     
    422436      ENDIF 
    423437 
    424       zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave 
    425438 
    426439      ! ----------------------------------------------------------------------------- ! 
     
    429442 
    430443      ! ... specific humidity at SST and IST tmask( 
    431       zsq(:,:) = 0.98 * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) ) 
     444      zsq(:,:) = rdct_qsat_salt * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) ) 
    432445      !! 
    433446      !! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate 
     
    444457      CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.5 
    445458         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    446       CASE( np_ECMWF     )   ;   CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! ECMWF 
    447          &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
     459 
     460      !LB: Skin!!! 
     461      CASE( np_ECMWF     ) 
     462         IF ( ln_skin ) THEN             
     463            !IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_ecmwf" WITH CSWL optional arrays!!!' 
     464            !IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => BEFORE ZST(40:50,30) =', zst(40:50,30) 
     465            CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! ECMWF          
     466               &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 
     467               &                                           Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1)) 
     468            !LB: "zst" and "zsq" have been updated with skin temp. !!! (from bulk SST)... 
     469            zst(:,:) = zst(:,:)*tmask(:,:,1) 
     470            zsq(:,:) = zsq(:,:)*tmask(:,:,1) 
     471            !IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => AFTER ZST(40:50,30) =', zst(40:50,30) 
     472         ELSE 
     473            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_ecmwf" WITHOUT CSWL optional arrays!!!' 
     474            CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! ECMWF          
     475               &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
     476         END IF 
     477         !LB. 
     478          
    448479      CASE DEFAULT 
    449480         CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 
    450481      END SELECT 
     482 
     483 
     484      !LB: cool-skin/warm-layer: 
     485      IF ( ln_skin ) tsk(:,:) = zst(:,:) 
     486 
    451487 
    452488      !                          ! Compute true air density : 
     
    507543 
    508544 
     545      ! ----------------------------------------------------------------------------- ! 
     546      !     III    Net longwave radiative FLUX                                        ! 
     547      ! ----------------------------------------------------------------------------- ! 
     548 
     549      !! LB: now after Turbulent fluxes because must use the skin temperature rather that the SST ! (zst is skin temperature if ln_skin==.TRUE.) 
     550      zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:,1) - emiss_w * stefan * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave 
     551 
     552 
    509553      IF(ln_ctl) THEN 
    510554         CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce: zqla   : ', tab2d_2=Ce_atm , clinfo2=' Ce_oce  : ' ) 
     
    519563 
    520564      ! ----------------------------------------------------------------------------- ! 
    521       !     III    Total FLUXES                                                       ! 
     565      !     IV    Total FLUXES                                                       ! 
    522566      ! ----------------------------------------------------------------------------- ! 
    523567      ! 
     
    527571      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar 
    528572         &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * rLfus                        &   ! remove latent melting heat for solid precip 
    529          &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST 
     573         &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST !LB??? pst is Celsius !? 
    530574         &     + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac  &   ! add liquid precip heat content at Tair 
    531575         &     * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp                          & 
     
    539583#endif 
    540584      ! 
    541       IF ( nn_ice == 0 ) THEN 
     585      !!#LB: NO WHY???? IF ( nn_ice == 0 ) THEN 
     586         CALL iom_put( "rho_air" ,   zrhoa )                ! output air density (kg/m^3) !#LB 
    542587         CALL iom_put( "qlw_oce" ,   zqlw )                 ! output downward longwave heat over the ocean 
    543588         CALL iom_put( "qsb_oce" , - zqsb )                 ! output downward sensible heat over the ocean 
     
    551596         CALL iom_put( 'snowpre', sprecip )                 ! Snow 
    552597         CALL iom_put( 'precip' , tprecip )                 ! Total precipitation 
    553       ENDIF 
     598         IF ( ln_skin ) THEN 
     599            CALL iom_put( "t_skin" ,  (zst - rt0) * tmask(:,:,1) )           ! T_skin in Celsius 
     600            CALL iom_put( "dt_skin" , (zst - pst - rt0) * tmask(:,:,1) )     ! T_skin - SST temperature difference... 
     601         END IF 
     602      !!#LB. ENDIF 
    554603      ! 
    555604      IF(ln_ctl) THEN 
     
    564613 
    565614 
    566  
    567    FUNCTION rho_air( ptak, pqa, pslp ) 
    568       !!------------------------------------------------------------------------------- 
    569       !!                           ***  FUNCTION rho_air  *** 
    570       !! 
    571       !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere 
    572       !! 
    573       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)  
    574       !!------------------------------------------------------------------------------- 
    575       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K] 
    576       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa       ! air specific humidity   [kg/kg] 
    577       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp      ! pressure in                [Pa] 
    578       REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air   ! density of moist air   [kg/m^3] 
    579       !!------------------------------------------------------------------------------- 
    580       ! 
    581       rho_air = pslp / (  R_dry*ptak * ( 1._wp + rctv0*pqa )  ) 
    582       ! 
    583    END FUNCTION rho_air 
    584  
    585  
    586    FUNCTION cp_air( pqa ) 
    587       !!------------------------------------------------------------------------------- 
    588       !!                           ***  FUNCTION cp_air  *** 
    589       !! 
    590       !! ** Purpose : Compute specific heat (Cp) of moist air 
    591       !! 
    592       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    593       !!------------------------------------------------------------------------------- 
    594       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg] 
    595       REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air   ! specific heat of moist air   [J/K/kg] 
    596       !!------------------------------------------------------------------------------- 
    597       ! 
    598       Cp_air = Cp_dry + Cp_vap * pqa 
    599       ! 
    600    END FUNCTION cp_air 
    601  
    602  
    603    FUNCTION q_sat( ptak, pslp ) 
    604       !!---------------------------------------------------------------------------------- 
    605       !!                           ***  FUNCTION q_sat  *** 
    606       !! 
    607       !! ** Purpose : Specific humidity at saturation in [kg/kg] 
    608       !!              Based on accurate estimate of "e_sat" 
    609       !!              aka saturation water vapor (Goff, 1957) 
    610       !! 
    611       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    612       !!---------------------------------------------------------------------------------- 
    613       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak    ! air temperature                       [K] 
    614       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp    ! sea level atmospheric pressure       [Pa] 
    615       REAL(wp), DIMENSION(jpi,jpj)             ::   q_sat   ! Specific humidity at saturation   [kg/kg] 
    616       ! 
    617       INTEGER  ::   ji, jj         ! dummy loop indices 
    618       REAL(wp) ::   ze_sat, ztmp   ! local scalar 
    619       !!---------------------------------------------------------------------------------- 
    620       ! 
    621       DO jj = 1, jpj 
    622          DO ji = 1, jpi 
    623             ! 
    624             ztmp = rt0 / ptak(ji,jj) 
    625             ! 
    626             ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957) 
    627             ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0)        & 
    628                &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) )  & 
    629                &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614  ) 
    630                ! 
    631             q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wp*pslp(ji,jj) - (1._wp - reps0)*ze_sat )   ! 0.01 because SLP is in [Pa] 
    632             ! 
    633          END DO 
    634       END DO 
    635       ! 
    636    END FUNCTION q_sat 
    637  
    638  
    639    FUNCTION gamma_moist( ptak, pqa ) 
    640       !!---------------------------------------------------------------------------------- 
    641       !!                           ***  FUNCTION gamma_moist  *** 
    642       !! 
    643       !! ** Purpose : Compute the moist adiabatic lapse-rate. 
    644       !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate 
    645       !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html 
    646       !! 
    647       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    648       !!---------------------------------------------------------------------------------- 
    649       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak          ! air temperature       [K] 
    650       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa           ! specific humidity [kg/kg] 
    651       REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist   ! moist adiabatic lapse-rate 
    652       ! 
    653       INTEGER  ::   ji, jj         ! dummy loop indices 
    654       REAL(wp) :: zrv, ziRT        ! local scalar 
    655       !!---------------------------------------------------------------------------------- 
    656       ! 
    657       DO jj = 1, jpj 
    658          DO ji = 1, jpi 
    659             zrv = pqa(ji,jj) / (1. - pqa(ji,jj)) 
    660             ziRT = 1. / (R_dry*ptak(ji,jj))    ! 1/RT 
    661             gamma_moist(ji,jj) = grav * ( 1. + rLevap*zrv*ziRT ) / ( Cp_dry + rLevap*rLevap*zrv*reps0*ziRT/ptak(ji,jj) ) 
    662          END DO 
    663       END DO 
    664       ! 
    665    END FUNCTION gamma_moist 
    666  
    667  
    668    FUNCTION L_vap( psst ) 
    669       !!--------------------------------------------------------------------------------- 
    670       !!                           ***  FUNCTION L_vap  *** 
    671       !! 
    672       !! ** Purpose : Compute the latent heat of vaporization of water from temperature 
    673       !! 
    674       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    675       !!---------------------------------------------------------------------------------- 
    676       REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap   ! latent heat of vaporization   [J/kg] 
    677       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K] 
    678       !!---------------------------------------------------------------------------------- 
    679       ! 
    680       L_vap = (  2.501 - 0.00237 * ( psst(:,:) - rt0)  ) * 1.e6 
    681       ! 
    682    END FUNCTION L_vap 
    683615 
    684616#if defined key_si3 
     
    710642      ! 
    711643      ! set transfer coefficients to default sea-ice values 
    712       Cd_atm(:,:) = Cd_ice 
    713       Ch_atm(:,:) = Cd_ice 
    714       Ce_atm(:,:) = Cd_ice 
     644      Cd_atm(:,:) = rCd_ice 
     645      Ch_atm(:,:) = rCd_ice 
     646      Ce_atm(:,:) = rCd_ice 
    715647 
    716648      wndm_ice(:,:) = 0._wp      !!gm brutal.... 
     
    737669      ENDIF 
    738670 
    739 !!      CALL iom_put( "Cd_ice", Cd_atm)  ! output value of pure ice-atm. transfer coef. 
     671!!      CALL iom_put( "rCd_ice", Cd_atm)  ! output value of pure ice-atm. transfer coef. 
    740672!!      CALL iom_put( "Ch_ice", Ch_atm)  ! output value of pure ice-atm. transfer coef. 
    741673 
     
    792724      REAL(wp) ::   zst3                     ! local variable 
    793725      REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      - 
    794       REAL(wp) ::   zztmp, z1_rLsub           !   -      - 
     726      REAL(wp) ::   zztmp, z1_rLsub          !   -      - 
    795727      REAL(wp) ::   zfr1, zfr2               ! local variables 
    796728      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_st         ! inverse of surface temperature 
     
    803735      !!--------------------------------------------------------------------- 
    804736      ! 
    805       zcoef_dqlw = 4.0 * 0.95 * Stef             ! local scalars 
    806       zcoef_dqla = -Ls * 11637800. * (-5897.8) 
     737      zcoef_dqlw = 4._wp * 0.95_wp * stefan             ! local scalars 
     738      zcoef_dqla = -rLsub * 11637800._wp * (-5897.8_wp) !LB: BAD! 
    807739      ! 
    808740      zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
     
    824756               qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj) 
    825757               ! Long  Wave (lw) 
    826                z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 
     758               z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - stefan * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 
    827759               ! lw sensitivity 
    828760               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3 
     
    834766               ! ... turbulent heat fluxes with Ch_atm recalculated in blk_ice_tau 
    835767               ! Sensible Heat 
    836                z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1)) 
     768               z_qsb(ji,jj,jl) = zrhoa(ji,jj) * rCp_air * Ch_atm(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1)) 
    837769               ! Latent Heat 
    838                qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls  * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  & 
     770               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * rLsub  * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  & 
    839771                  &                ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) ) 
    840772               ! Latent heat sensitivity for ice (Dqla/Dt) 
     
    847779 
    848780               ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) 
    849                z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj) 
     781               z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * rCp_air * Ch_atm(ji,jj) * wndm_ice(ji,jj) 
    850782 
    851783               ! ----------------------------! 
     
    1064996      ! generic drag over a cell partly covered by ice 
    1065997      !!Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) +  &                        ! pure ocean drag 
    1066       !!   &      Cd_ice      *           at_i_b(:,:)   +  &                        ! pure ice drag 
     998      !!   &      rCd_ice      *           at_i_b(:,:)   +  &                        ! pure ice drag 
    1067999      !!   &      zCe         * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**zmu   ! change due to sea-ice morphology 
    10681000 
    10691001      ! ice-atm drag 
    1070       Cd(:,:) = Cd_ice +  &                                                         ! pure ice drag 
     1002      Cd(:,:) = rCd_ice +  &                                                         ! pure ice drag 
    10711003         &      zCe    * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp)  ! change due to sea-ice morphology 
    10721004       
     
    11411073      ! Atmospheric and Surface Variables 
    11421074      zst(:,:)     = sst_m(:,:) + rt0                                        ! convert SST from Celcius to Kelvin 
    1143       zqo_sat(:,:) = 0.98_wp * q_sat( zst(:,:)   , sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ocean [kg/kg] 
    1144       zqi_sat(:,:) = 0.98_wp * q_sat( ztm_su(:,:), sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ice   [kg/kg] 
     1075      zqo_sat(:,:) = rdct_qsat_salt * q_sat( zst(:,:)   , sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ocean [kg/kg] 
     1076      zqi_sat(:,:) =                  q_sat( ztm_su(:,:), sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ice   [kg/kg] !LB: no 0.98 !!(rdct_qsat_salt) 
    11451077      ! 
    11461078      DO jj = 2, jpjm1           ! reduced loop is necessary for reproducibility 
Note: See TracChangeset for help on using the changeset viewer.