Changeset 11111


Ignore:
Timestamp:
2019-06-14T15:55:28+02:00 (17 months 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".

Location:
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE
Files:
2 added
7 edited

Legend:

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

    r10068 r11111  
    5050   REAL(wp), PUBLIC ::   sice     =    6.0_wp        !: salinity of ice (for pisces)          [psu] 
    5151   REAL(wp), PUBLIC ::   soce     =   34.7_wp        !: salinity of sea (for pisces and isf)  [psu] 
    52    REAL(wp), PUBLIC ::   rLevap   =    2.5e+6_wp     !: latent heat of evaporation (water) 
     52   REAL(wp), PUBLIC ::   rLevap   =    2.46e+6_wp    !: latent heat of vaporization for sea-water   [J/kg] #LB 
    5353   REAL(wp), PUBLIC ::   vkarmn   =    0.4_wp        !: von Karman constant 
    5454   REAL(wp), PUBLIC ::   stefan   =    5.67e-8_wp    !: Stefan-Boltzmann constant  
     
    6666   REAL(wp), PUBLIC ::   r1_rhos                     !: 1 / rhos 
    6767   REAL(wp), PUBLIC ::   r1_rcpi                     !: 1 / rcpi 
     68 
     69 
     70   !#LB: 
     71   !! Mainly used in OCE/SBC (removed from sbcblk.F90) 
     72   REAL(wp), PARAMETER, PUBLIC :: rCp_dry = 1005.0_wp   !: Specic heat of dry air, constant pressure      [J/K/kg] 
     73   REAL(wp), PARAMETER, PUBLIC :: rCp_vap = 1860.0_wp   !: Specic heat of water vapor, constant pressure  [J/K/kg] 
     74   REAL(wp), PARAMETER, PUBLIC :: R_dry   = 287.05_wp   !: Specific gas constant for dry air              [J/K/kg] 
     75   REAL(wp), PARAMETER, PUBLIC :: R_vap   = 461.495_wp  !: Specific gas constant for water vapor          [J/K/kg] 
     76   REAL(wp), PARAMETER, PUBLIC :: reps0   = R_dry/R_vap !: ratio of gas constant for dry air and water vapor => ~ 0.622 
     77   REAL(wp), PARAMETER, PUBLIC :: rctv0   = R_vap/R_dry !: for virtual temperature (== (1-eps)/eps) => ~ 0.608 
     78   REAL(wp), PARAMETER, PUBLIC :: rCp_air = 1000.5_wp   !: specific heat of air (only used for ice fluxes now...) 
     79   REAL(wp), PARAMETER, PUBLIC :: rCd_ice = 1.4e-3_wp   !: transfer coefficient over ice 
     80   REAL(wp), PARAMETER, PUBLIC :: albo    = 0.066_wp    !: ocean albedo assumed to be constant 
     81   ! 
     82   REAL(wp), PARAMETER, PUBLIC :: rho0_a  = 1.2_wp      !: Approx. of density of air                       [kg/m^3] 
     83   REAL(wp), PARAMETER, PUBLIC :: rho0_w  = 1025._wp    !: Density of sea-water  (ECMWF->1025)             [kg/m^3] 
     84   REAL(wp), PARAMETER, PUBLIC :: rCp0_w  = 4190._wp    !: Specific heat capacity of seawater (ECMWF 4190) [J/K/kg] 
     85   REAL(wp), PARAMETER, PUBLIC :: rnu0_w  = 1.e-6_wp    !: kinetic viscosity of water                      [m^2/s] 
     86   REAL(wp), PARAMETER, PUBLIC :: rk0_w   = 0.6_wp      !: thermal conductivity of water (at 20C)          [W/m/K] 
     87   ! 
     88   REAL(wp), PARAMETER, PUBLIC :: emiss_w = 1._wp       !: Surface emissivity (black-body long-wave radiation) of sea-water [] 
     89   !                                                    !: Theoretically close to 0.97! Yet, taken equal as 1 to account for 
     90   !                                                    !: the small fraction of downwelling shortwave reflected at the 
     91   !                                                    !: surface (Lind & Katsaros, 1986) 
     92   REAL(wp), PARAMETER, PUBLIC :: rdct_qsat_salt = 0.98_wp  !: reduction factor on specific humidity at saturation (q_sat(T_s)) due to salt   
     93   !#LB. 
     94 
     95 
     96    
    6897   !!---------------------------------------------------------------------- 
    6998   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbc_oce.F90

    r10882 r11111  
    149149   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   frq_m     !: mean (nn_fsbc time-step) fraction of solar net radiation absorbed in the 1st T level [-] 
    150150 
     151   !!---------------------------------------------------------------------- 
     152   !!                     Cool-skin/Warm-layer 
     153   !!---------------------------------------------------------------------- 
     154   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tsk       !: sea-surface skin temperature out of the cool-skin/warm-layer parameterization [Celsius] 
     155 
     156    
    151157   !! * Substitutions 
    152158#  include "vectopt_loop_substitute.h90" 
  • 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 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare.F90

    r10069 r11111  
    11MODULE sbcblk_algo_coare 
    22   !!====================================================================== 
    3    !!                       ***  MODULE  sbcblk_algo_coare  *** 
    4    !! Computes turbulent components of surface fluxes 
    5    !!         according to Fairall et al. 2003 (COARE v3) 
    6    !! 
     3   !!                   ***  MODULE  sbcblk_algo_coare  *** 
     4   !! Computes: 
    75   !!   * bulk transfer coefficients C_D, C_E and C_H 
    86   !!   * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed 
     
    108   !!   => all these are used in bulk formulas in sbcblk.F90 
    119   !! 
    12    !!    Using the bulk formulation/param. of COARE v3, Fairall et al. 2003 
    13    !! 
     10   !!    Using the bulk formulation/param. of COARE v3, Fairall et al., 2003 
    1411   !! 
    1512   !!       Routine turb_coare maintained and developed in AeroBulk 
    16    !!                     (http://aerobulk.sourceforge.net/) 
     13   !!                     (https://github.com/brodeau/aerobulk) 
    1714   !! 
    18    !!            Author: Laurent Brodeau, 2016, brodeau@gmail.com 
    19    !! 
    20    !!====================================================================== 
    21    !! History :  3.6  !  2016-02  (L.Brodeau)   Original code 
     15   !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk) 
     16   !!---------------------------------------------------------------------- 
     17   !! History :  4.0  !  2016-02  (L.Brodeau)   Original code 
    2218   !!---------------------------------------------------------------------- 
    2319 
     
    2723   !!                   returns the effective bulk wind speed at 10m 
    2824   !!---------------------------------------------------------------------- 
    29    USE oce            ! ocean dynamics and tracers 
    30    USE dom_oce        ! ocean space and time domain 
    31    USE phycst         ! physical constants 
    32    USE sbc_oce        ! Surface boundary condition: ocean fields 
    33    USE sbcwave, ONLY  :  cdn_wave ! wave module 
     25   USE oce             ! ocean dynamics and tracers 
     26   USE dom_oce         ! ocean space and time domain 
     27   USE phycst          ! physical constants 
     28   USE sbc_oce         ! Surface boundary condition: ocean fields 
     29   USE sbcwave, ONLY   :  cdn_wave ! wave module 
    3430#if defined key_si3 || defined key_cice 
    35    USE sbc_ice        ! Surface boundary condition: ice fields 
     31   USE sbc_ice         ! Surface boundary condition: ice fields 
    3632#endif 
    37    ! 
     33   USE sbcblk_phymbl  !LB: all thermodynamics functions, rho_air, q_sat, etc... 
    3834   USE in_out_manager ! I/O manager 
    3935   USE iom            ! I/O manager library 
     
    5046   REAL(wp), PARAMETER ::   zi0     = 600._wp      ! scale height of the atmospheric boundary layer... 
    5147   REAL(wp), PARAMETER ::   Beta0   =   1.250_wp   ! gustiness parameter 
    52    REAL(wp), PARAMETER ::   rctv0   =   0.608_wp   ! constant to obtain virtual temperature... 
     48 
     49   INTEGER , PARAMETER ::   nb_itt = 5             ! number of itterations 
    5350 
    5451   !!---------------------------------------------------------------------- 
     
    6158      !!                      ***  ROUTINE  turb_coare  *** 
    6259      !! 
    63       !!            2015: L. Brodeau (brodeau@gmail.com) 
    64       !! 
    6560      !! ** Purpose :   Computes turbulent transfert coefficients of surface 
    6661      !!                fluxes according to Fairall et al. (2003) 
    6762      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu 
    68       !! 
    69       !! ** Method : Monin Obukhov Similarity Theory 
    70       !!---------------------------------------------------------------------- 
     63      !!                Returns the effective bulk wind speed at 10m to be used in the bulk formulas 
     64      !! 
    7165      !! 
    7266      !! INPUT : 
     
    8882      !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K] 
    8983      !!    *  q_zu   : specific humidity of air        //                    [kg/kg] 
    90       !!    *  U_blk  : bulk wind at 10m                                      [m/s] 
    91       !!---------------------------------------------------------------------- 
     84      !!    *  U_blk  : bulk wind speed at 10m                                [m/s] 
     85      !! 
     86      !! 
     87      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk) 
     88      !!---------------------------------------------------------------------------------- 
    9289      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m] 
    9390      REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                             [m] 
     
    106103      ! 
    107104      INTEGER :: j_itt 
    108       LOGICAL ::   l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U 
    109       INTEGER , PARAMETER ::   nb_itt = 4       ! number of itterations 
     105      LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U 
    110106 
    111107      REAL(wp), DIMENSION(jpi,jpj) ::  & 
     
    125121 
    126122      !! First guess of temperature and humidity at height zu: 
    127       t_zu = MAX(t_zt , 0.0)    ! who knows what's given on masked-continental regions... 
    128       q_zu = MAX(q_zt , 1.e-6)  !               " 
     123      t_zu = MAX( t_zt , 199.0_wp )   ! who knows what's given on masked-continental regions... 
     124      q_zu = MAX( q_zt , 1.e-6_wp )   !               " 
    129125 
    130126      !! Pot. temp. difference (and we don't want it to be 0!) 
    131       dt_zu = t_zu - sst ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu ) 
    132       dq_zu = q_zu - ssq ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu ) 
    133  
    134       znu_a = visc_air(t_zt) ! Air viscosity (m^2/s) at zt given from temperature in (K) 
     127      dt_zu = t_zu - sst ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 
     128      dq_zu = q_zu - ssq ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 
     129 
     130      znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K) 
    135131 
    136132      ztmp2 = 0.5*0.5  ! initial guess for wind gustiness contribution 
     
    144140 
    145141      z0     = alfa_charn(U_blk)*u_star*u_star/grav + 0.11*znu_a/u_star 
    146       z0t    = 0.1*EXP(vkarmn/(0.00115/(vkarmn/ztmp1)))   !  WARNING: 1/z0t ! 
     142      z0     = MIN(ABS(z0), 0.001)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO 
     143      z0t    = 1. / ( 0.1*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) ) 
     144      z0t    = MIN(ABS(z0t), 0.001)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO 
    147145 
    148146      ztmp2  = vkarmn/ztmp0 
    149147      Cd     = ztmp2*ztmp2    ! first guess of Cd 
    150148 
    151       ztmp0 = vkarmn*vkarmn/LOG(zt*z0t)/Cd 
    152  
    153       !Ribcu = -zu/(zi0*0.004*Beta0**3) !! Saturation Rib, zi0 = tropicalbound. layer depth 
    154       ztmp2  = grav*zu*(dt_zu + rctv0*t_zu*dq_zu)/(t_zu*U_blk*U_blk)  !! Ribu Bulk Richardson number 
    155       ztmp1 = 0.5 + sign(0.5 , ztmp2) 
     149      ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd 
     150 
     151      ztmp2  = grav*zu*(dt_zu + rctv0*t_zu*dq_zu)/(t_zu*U_blk*U_blk)  !! Ribu Bulk Richardson number ;       !Ribcu = -zu/(zi0*0.004*Beta0**3) !! Saturation Rib, zi0 = tropicalbound. layer depth 
     152 
     153      !! First estimate of zeta_u, depending on the stability, ie sign of Ribu (ztmp2): 
     154      ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 ) 
    156155      ztmp0 = ztmp0*ztmp2 
    157       !!             Ribu < 0                                 Ribu > 0   Beta = 1.25 
    158       zeta_u = (1.-ztmp1) * (ztmp0/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) & 
    159          &  +     ztmp1   * (ztmp0*(1. + 27./9.*ztmp2/ztmp0)) 
     156      zeta_u = (1.-ztmp1) * (ztmp0/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) & !  Ribu < 0 
     157         &  +     ztmp1   * (ztmp0*(1. + 27./9.*ztmp2/ztmp0))               !  Ribu > 0 
     158      !#LOLO: should make sure that the "ztmp0" of "27./9.*ztmp2/ztmp0" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" ! 
    160159 
    161160      !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L 
    162       ztmp0  =  vkarmn/(LOG(zu*z0t) - psi_h_coare(zeta_u)) 
     161      ztmp0  = vkarmn/(LOG(zu/z0t) - psi_h_coare(zeta_u)) 
    163162 
    164163      u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_coare(zeta_u)) 
     
    168167      ! What's need to be done if zt /= zu: 
    169168      IF( .NOT. l_zt_equal_zu ) THEN 
    170  
     169         !! First update of values at zu (or zt for wind) 
    171170         zeta_t = zt*zeta_u/zu 
    172  
    173          !! First update of values at zu (or zt for wind) 
    174171         ztmp0 = psi_h_coare(zeta_u) - psi_h_coare(zeta_t) 
    175          ztmp1 = log(zt/zu) + ztmp0 
     172         ztmp1 = LOG(zt/zu) + ztmp0 
    176173         t_zu = t_zt - t_star/vkarmn*ztmp1 
    177174         q_zu = q_zt - q_star/vkarmn*ztmp1 
    178          q_zu = (0.5 + sign(0.5,q_zu))*q_zu !Makes it impossible to have negative humidity : 
    179  
    180          dt_zu = t_zu - sst  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu ) 
    181          dq_zu = q_zu - ssq  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu ) 
    182  
     175         q_zu = (0.5 + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity : 
     176         ! 
     177         dt_zu = t_zu - sst  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 
     178         dq_zu = q_zu - ssq  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 
    183179      END IF 
    184180 
     
    188184         !!Inverse of Monin-Obukov length (1/L) : 
    189185         ztmp0 = One_on_L(t_zu, q_zu, u_star, t_star, q_star)  ! 1/L == 1/[Monin-Obukhov length] 
     186         ztmp0 = SIGN( MIN(ABS(ztmp0),200._wp), ztmp0 ) ! (prevents FPE from stupid values from masked region later on...) !#LOLO 
    190187 
    191188         ztmp1 = u_star*u_star   ! u*^2 
     
    193190         !! Update wind at 10m taking into acount convection-related wind gustiness: 
    194191         ! Ug = Beta*w*  (Beta = 1.25, Fairall et al. 2003, Eq.8): 
    195          ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0.))**(2./3.)   ! => ztmp2 == Ug^2 
     192         ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0._wp))**(2./3.)   ! => ztmp2 == Ug^2 
    196193         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before 600. 
    197          U_blk = MAX(sqrt(U_zu*U_zu + ztmp2), 0.2)        ! include gustiness in bulk wind speed 
     194         U_blk = MAX(sqrt(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed 
    198195         ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0. 
    199196 
     
    203200         !! Roughness lengthes z0, z0t (z0q = z0t) : 
    204201         z0   = ztmp2*ztmp1/grav + 0.11*znu_a/u_star ! Roughness length (eq.6) 
    205          ztmp1 = z0*u_star/znu_a                             ! Re_r: roughness Reynolds number 
    206          z0t  = min( 1.1E-4 , 5.5E-5*ztmp1**(-0.6) )        ! Scalar roughness for both theta and q (eq.28) 
     202         ztmp1 = z0*u_star/znu_a                           ! Re_r: roughness Reynolds number 
     203         z0t  = MIN( 1.1E-4_wp , 5.5E-5_wp*ztmp1**(-0.6_wp) ) ! Scalar roughness for both theta and q (eq.28) 
    207204 
    208205         !! Stability parameters: 
    209          zeta_u = zu*ztmp0 ; zeta_u = sign( min(abs(zeta_u),50.0), zeta_u ) 
     206         zeta_u = zu*ztmp0 
     207         zeta_u = SIGN( MIN(ABS(zeta_u),50.0_wp), zeta_u ) 
    210208         IF( .NOT. l_zt_equal_zu ) THEN 
    211             zeta_t = zt*ztmp0 ;  zeta_t = sign( min(abs(zeta_t),50.0), zeta_t ) 
     209            zeta_t = zt*ztmp0 
     210            zeta_t = SIGN( MIN(ABS(zeta_t),50.0_wp), zeta_t ) 
    212211         END IF 
    213212 
    214213         !! Turbulent scales at zu=10m : 
    215214         ztmp0   = psi_h_coare(zeta_u) 
    216          ztmp1   = vkarmn/(LOG(zu) -LOG(z0t) - ztmp0) 
     215         ztmp1   = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0) 
    217216 
    218217         t_star = dt_zu*ztmp1 
    219218         q_star = dq_zu*ztmp1 
    220          u_star = U_blk*vkarmn/(LOG(zu) -LOG(z0) - psi_m_coare(zeta_u)) 
     219         u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_coare(zeta_u)) 
    221220 
    222221         IF( .NOT. l_zt_equal_zu ) THEN 
     
    259258      !! Wind greater than 18 m/s :  alfa = 0.018 
    260259      !! 
    261       !! Author: L. Brodeau, june 2016 / AeroBulk  (https://sourceforge.net/p/aerobulk) 
     260      !! Author: L. Brodeau, june 2016 / AeroBulk  (https://github.com/brodeau/aerobulk/) 
    262261      !!------------------------------------------------------------------- 
    263262      REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn 
     
    274273            ! 
    275274            ! Charnock's constant, increases with the wind : 
    276             zgt10 = 0.5 + SIGN(0.5,(zw - 10.)) ! If zw<10. --> 0, else --> 1 
    277             zgt18 = 0.5 + SIGN(0.5,(zw - 18.)) ! If zw<18. --> 0, else --> 1 
     275            zgt10 = 0.5 + SIGN(0.5_wp,(zw - 10)) ! If zw<10. --> 0, else --> 1 
     276            zgt18 = 0.5 + SIGN(0.5_wp,(zw - 18.)) ! If zw<18. --> 0, else --> 1 
    278277            ! 
    279278            alfa_charn(ji,jj) =  (1. - zgt10)*0.011    &    ! wind is lower than 10 m/s 
     
    294293      !! 
    295294      !! Author: L. Brodeau, june 2016 / AeroBulk 
    296       !!         (https://sourceforge.net/p/aerobulk) 
     295      !!         (https://github.com/brodeau/aerobulk/) 
    297296      !!------------------------------------------------------------------------ 
    298297      REAL(wp), DIMENSION(jpi,jpj)             :: One_on_L         !: 1./(Monin Obukhov length) [m^-1] 
     
    330329      !!       form from Beljaars and Holtslag (1991) 
    331330      !! 
    332       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     331      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    333332      !!---------------------------------------------------------------------------------- 
    334333      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare 
     
    356355            zf = zta*zta 
    357356            zf = zf/(1. + zf) 
    358             zc = MIN(50., 0.35*zta) 
    359             zstab = 0.5 + SIGN(0.5, zta) 
     357            zc = MIN(50._wp, 0.35_wp*zta) 
     358            zstab = 0.5 + SIGN(0.5_wp, zta) 
    360359            ! 
    361360            psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0) 
     
    383382      !! 
    384383      !! Author: L. Brodeau, june 2016 / AeroBulk 
    385       !!         (https://sourceforge.net/p/aerobulk) 
     384      !!         (https://github.com/brodeau/aerobulk/) 
    386385      !!---------------------------------------------------------------- 
    387386      !! 
     
    408407            zf = zta*zta 
    409408            zf = zf/(1. + zf) 
    410             zc = MIN(50.,0.35*zta) 
    411             zstab = 0.5 + SIGN(0.5, zta) 
     409            zc = MIN(50._wp,0.35_wp*zta) 
     410            zstab = 0.5 + SIGN(0.5_wp, zta) 
    412411            ! 
    413412            psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & 
     
    420419   END FUNCTION psi_h_coare 
    421420 
    422  
    423    FUNCTION visc_air( ptak ) 
    424       !!--------------------------------------------------------------------- 
    425       !! Air kinetic viscosity (m^2/s) given from temperature in degrees... 
    426       !! 
    427       !! Author: L. Brodeau, june 2016 / AeroBulk 
    428       !!         (https://sourceforge.net/p/aerobulk) 
    429       !!--------------------------------------------------------------------- 
    430       !! 
    431       REAL(wp), DIMENSION(jpi,jpj) :: visc_air 
    432       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak  ! air temperature in (K) 
    433       ! 
    434       INTEGER  ::   ji, jj         ! dummy loop indices 
    435       REAL(wp) ::   ztc, ztc2      ! local scalar 
    436       ! 
    437       DO jj = 1, jpj 
    438          DO ji = 1, jpi 
    439             ztc  = ptak(ji,jj) - rt0   ! air temp, in deg. C 
    440             ztc2 = ztc*ztc 
    441             visc_air(ji,jj) = 1.326E-5*(1. + 6.542E-3*ztc + 8.301E-6*ztc2 - 4.84E-9*ztc2*ztc) 
    442          END DO 
    443       END DO 
    444       ! 
    445    END FUNCTION visc_air 
    446  
    447  
     421   !!====================================================================== 
    448422END MODULE sbcblk_algo_coare 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p5.F90

    r10069 r11111  
    1414   !! 
    1515   !!       Routine turb_coare3p5 maintained and developed in AeroBulk 
    16    !!                     (http://aerobulk.sourceforge.net/) 
     16   !!                     (https://github.com/brodeau/aerobulk/) 
    1717   !! 
    1818   !!            Author: Laurent Brodeau, 2016, brodeau@gmail.com 
     
    5151   REAL(wp), PARAMETER ::   zi0        = 600.      ! scale height of the atmospheric boundary layer...1 
    5252   REAL(wp), PARAMETER ::   Beta0      =   1.25    ! gustiness parameter 
    53    REAL(wp), PARAMETER ::   rctv0      =   0.608   ! constant to obtain virtual temperature... 
    5453 
    5554   !!---------------------------------------------------------------------- 
     
    6867      !! ** Method : Monin Obukhov Similarity Theory 
    6968      !! 
    70       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)  
     69      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)  
    7170      !! 
    7271      !! INPUT : 
     
    265264      !! 
    266265      !! Author: L. Brodeau, june 2016 / AeroBulk 
    267       !!         (https://sourceforge.net/p/aerobulk) 
     266      !!         (https://github.com/brodeau/aerobulk/) 
    268267      !!------------------------------------------------------------------------ 
    269268      REAL(wp), DIMENSION(jpi,jpj)             :: One_on_L         !: 1./(Monin Obukhov length) [m^-1] 
     
    300299      !!       form from Beljaars and Holtslag (1991) 
    301300      !! 
    302       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     301      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    303302      !!---------------------------------------------------------------------------------- 
    304303      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare 
     
    353352      !! 
    354353      !! Author: L. Brodeau, june 2016 / AeroBulk 
    355       !!         (https://sourceforge.net/p/aerobulk) 
     354      !!         (https://github.com/brodeau/aerobulk/) 
    356355      !!---------------------------------------------------------------- 
    357356      !! 
     
    396395      !! 
    397396      !! Author: L. Brodeau, june 2016 / AeroBulk 
    398       !!         (https://sourceforge.net/p/aerobulk) 
     397      !!         (https://github.com/brodeau/aerobulk/) 
    399398      !!--------------------------------------------------------------------- 
    400399      REAL(wp), DIMENSION(jpi,jpj)             ::   visc_air 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90

    r10069 r11111  
    11MODULE sbcblk_algo_ecmwf 
    22   !!====================================================================== 
    3    !!                       ***  MODULE  sbcblk_algo_ecmwf  *** 
    4    !! Computes turbulent components of surface fluxes 
    5    !!         according to the method in IFS of the ECMWF model 
    6    !! 
     3   !!                   ***  MODULE  sbcblk_algo_ecmwf  *** 
     4   !! Computes: 
    75   !!   * bulk transfer coefficients C_D, C_E and C_H 
    86   !!   * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed 
     
    108   !!   => all these are used in bulk formulas in sbcblk.F90 
    119   !! 
    12    !!    Using the bulk formulation/param. of IFS of ECMWF (cycle 31r2) 
     10   !!    Using the bulk formulation/param. of IFS of ECMWF (cycle 40r1) 
    1311   !!         based on IFS doc (avaible online on the ECMWF's website) 
    1412   !! 
     13   !!       Routine turb_ecmwf maintained and developed in AeroBulk 
     14   !!                     (https://github.com/brodeau/aerobulk) 
    1515   !! 
    16    !!       Routine turb_ecmwf maintained and developed in AeroBulk 
    17    !!                     (http://aerobulk.sourceforge.net/) 
    18    !! 
    19    !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     16   !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk) 
    2017   !!---------------------------------------------------------------------- 
    2118   !! History :  4.0  !  2016-02  (L.Brodeau)   Original code 
     
    4138 
    4239   USE sbc_oce         ! Surface boundary condition: ocean fields 
     40   USE sbcblk_phymbl   !LB: all thermodynamics functions, rho_air, q_sat, etc... 
     41   USE sbcblk_skin     ! cool-skin/warm layer scheme (CSWL_ECMWF) 
    4342 
    4443   IMPLICIT NONE 
     
    5251   REAL(wp), PARAMETER ::   zi0     = 1000.   ! scale height of the atmospheric boundary layer...1 
    5352   REAL(wp), PARAMETER ::   Beta0    = 1.     ! gustiness parameter ( = 1.25 in COAREv3) 
    54    REAL(wp), PARAMETER ::   rctv0    = 0.608  ! constant to obtain virtual temperature... 
    55    REAL(wp), PARAMETER ::   Cp_dry = 1005.0   ! Specic heat of dry air, constant pressure      [J/K/kg] 
    56    REAL(wp), PARAMETER ::   Cp_vap = 1860.0   ! Specic heat of water vapor, constant pressure  [J/K/kg] 
    5753   REAL(wp), PARAMETER ::   alpha_M = 0.11    ! For roughness length (smooth surface term) 
    5854   REAL(wp), PARAMETER ::   alpha_H = 0.40    ! (Chapter 3, p.34, IFS doc Cy31r1) 
    5955   REAL(wp), PARAMETER ::   alpha_Q = 0.62    ! 
     56 
     57   INTEGER , PARAMETER ::   nb_itt = 5        ! number of itterations 
     58 
     59   !! Cool-Skin / Warm-Layer related parameters: 
     60   REAL(wp),    PARAMETER :: & 
     61      &  rdt0    = 3600.*1.5, & !: time step 
     62      &  rd0     = 3. ,      &  !: Depth scale [m], "d" in Eq.11 (Zeng & Beljaars 2005) 
     63      &  rNu0    = 0.5          !: Nu (exponent of temperature profile) Eq.11 
     64   !                            !: (Zeng & Beljaars 2005) !: set to 0.5 instead of 
     65   !                            !: 0.3 to respect a warming of +3 K in calm 
     66   !                            !: condition for the insolation peak of +1000W/m^2 
     67   INTEGER,    PARAMETER :: & 
     68      &  nb_itt_wl = 10         !: number of sub-itterations for solving the differential equation in warm-layer part 
     69   !                            !:  => use "nb_itt_wl = 1" for No itteration! => way cheaper !!! 
     70   !                            !:    => assumes balance between the last 2 terms of Eq.11 (lhs of eq.11 = 0) 
     71   !                            !:    => in that case no need for sub-itterations ! 
     72   !                            !:    => ACCEPTABLE IN MOST CONDITIONS ! (UNLESS: sunny + very calm/low-wind conditions) 
     73   !                            !:  => Otherwize use "nb_itt_wl = 10" 
     74 
     75 
     76 
     77 
    6078   !!---------------------------------------------------------------------- 
    6179CONTAINS 
    6280 
    63    SUBROUTINE TURB_ECMWF( zt, zu, sst, t_zt, ssq , q_zt , U_zu,   & 
    64       &                   Cd, Ch, Ce , t_zu, q_zu, U_blk,         & 
    65       &                   Cdn, Chn, Cen                           ) 
     81   SUBROUTINE TURB_ECMWF( zt, zu, T_s, t_zt, q_s, q_zt, U_zu, & 
     82      &                   Cd, Ch, Ce, t_zu, q_zu, U_blk,      & 
     83      &                   Cdn, Chn, Cen,                      & 
     84      &                   Qsw, rad_lw, slp                   ) 
    6685      !!---------------------------------------------------------------------------------- 
    6786      !!                      ***  ROUTINE  turb_ecmwf  *** 
    6887      !! 
    69       !!            2015: L. Brodeau (brodeau@gmail.com) 
    70       !! 
    7188      !! ** Purpose :   Computes turbulent transfert coefficients of surface 
    72       !!                fluxes according to IFS doc. (cycle 31) 
     89      !!                fluxes according to IFS doc. (cycle 40) 
    7390      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu 
    74       !! 
    75       !! ** Method : Monin Obukhov Similarity Theory 
     91      !!                Returns the effective bulk wind speed at 10m to be used in the bulk formulas 
     92      !! 
     93      !!                Applies the cool-skin warm-layer correction of the SST to T_s 
     94      !!                if the net shortwave flux at the surface (Qsw), the downwelling longwave 
     95      !!                radiative fluxes at the surface (rad_lw), and the sea-leve pressure (slp) 
     96      !!                are provided as (optional) arguments! 
    7697      !! 
    7798      !! INPUT : 
     
    80101      !!    *  zu   : height for wind speed (generally 10m)                   [m] 
    81102      !!    *  U_zu : scalar wind speed at 10m                                [m/s] 
    82       !!    *  sst  : SST                                                     [K] 
    83103      !!    *  t_zt : potential air temperature at zt                         [K] 
    84       !!    *  ssq  : specific humidity at saturation at SST                  [kg/kg] 
    85104      !!    *  q_zt : specific humidity of air at zt                          [kg/kg] 
    86105      !! 
     106      !! INPUT/OUTPUT: 
     107      !! ------------- 
     108      !!    *  T_s  : SST or skin temperature                                 [K] 
     109      !!    *  q_s  : SSQ aka saturation specific humidity at temp. T_s       [kg/kg] 
     110      !!              -> doesn't need to be given a value if skin temp computed (in case l_use_skin=True) 
     111      !!              -> MUST be given the correct value if not computing skint temp. (in case l_use_skin=False) 
     112      !! 
     113      !! OPTIONAL INPUT (will trigger l_use_skin=TRUE if present!): 
     114      !! --------------- 
     115      !!    *  Qsw    : net solar flux (after albedo) at the surface (>0)     [W/m^2] 
     116      !!    *  rad_lw : downwelling longwave radiation at the surface  (>0)   [W/m^2] 
     117      !!    *  slp    : sea-level pressure                                    [Pa] 
    87118      !! 
    88119      !! OUTPUT : 
     
    93124      !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K] 
    94125      !!    *  q_zu   : specific humidity of air        //                    [kg/kg] 
    95       !!    *  U_blk  : bulk wind at 10m                                      [m/s] 
    96       !! 
    97       !! 
    98       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     126      !!    *  U_blk  : bulk wind speed at 10m                                [m/s] 
     127      !! 
     128      !! 
     129      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    99130      !!---------------------------------------------------------------------------------- 
    100131      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m] 
    101132      REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                             [m] 
    102       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   sst      ! sea surface temperature                [Kelvin] 
     133      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) ::   T_s      ! sea surface temperature                [Kelvin] 
    103134      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   t_zt     ! potential air temperature              [Kelvin] 
    104       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   ssq      ! sea surface specific humidity           [kg/kg] 
    105       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity                   [kg/kg] 
     135      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) ::   q_s      ! sea surface specific humidity           [kg/kg] 
     136      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity at zt             [kg/kg] 
    106137      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu                [m/s] 
    107138      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau) 
     
    113144      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients 
    114145      ! 
     146      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   Qsw      !             [W/m^2] 
     147      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   rad_lw   !             [W/m^2] 
     148      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   slp      !             [Pa] 
     149      ! 
    115150      INTEGER :: j_itt 
    116       LOGICAL ::   l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U 
    117       INTEGER , PARAMETER ::   nb_itt = 4       ! number of itterations 
    118       ! 
    119       REAL(wp), DIMENSION(jpi,jpj) ::   u_star, t_star, q_star,  & 
     151      LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U 
     152      ! 
     153      REAL(wp), DIMENSION(jpi,jpj) ::  & 
     154         &  u_star, t_star, q_star, & 
    120155         &  dt_zu, dq_zu,    & 
    121156         &  znu_a,           & !: Nu_air, Viscosity of air 
    122157         &  Linv,            & !: 1/L (inverse of Monin Obukhov length... 
    123158         &  z0, z0t, z0q 
     159      ! 
     160      ! Cool skin: 
     161      LOGICAL :: l_use_skin = .FALSE. 
     162      REAL(wp), DIMENSION(jpi,jpj) :: zsst    ! to back up the initial bulk SST 
     163 
     164      ! 
    124165      REAL(wp), DIMENSION(jpi,jpj) ::   func_m, func_h 
    125166      REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2 
    126167      !!---------------------------------------------------------------------------------- 
    127168      ! 
     169      ! Cool skin ? 
     170      IF( PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) ) THEN 
     171         l_use_skin = .TRUE. 
     172      END IF 
     173      IF (lwp) PRINT *, ' *** LOLO(sbcblk_algo_ecmwf.F90) => l_use_skin =', l_use_skin 
     174 
    128175      ! Identical first gess as in COARE, with IFS parameter values though 
    129176      ! 
     
    131178      IF( ABS(zu - zt) < 0.01 )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision 
    132179 
     180      !! Initialization for cool skin: 
     181      IF( l_use_skin ) THEN 
     182         zsst   = T_s    ! save the bulk SST 
     183         T_s    = T_s - 0.25                      ! First guess of correction 
     184         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s 
     185      END IF 
    133186 
    134187      !! First guess of temperature and humidity at height zu: 
    135       t_zu = MAX( t_zt , 0.0  )   ! who knows what's given on masked-continental regions... 
    136       q_zu = MAX( q_zt , 1.e-6)   !               " 
     188      t_zu = MAX( t_zt , 0.0_wp  )   ! who knows what's given on masked-continental regions... 
     189      q_zu = MAX( q_zt , 1.e-6_wp)   !               " 
    137190 
    138191      !! Pot. temp. difference (and we don't want it to be 0!) 
    139       dt_zu = t_zu - sst   ;   dt_zu = SIGN( MAX(ABS(dt_zu),1.e-6), dt_zu ) 
    140       dq_zu = q_zu - ssq   ;   dq_zu = SIGN( MAX(ABS(dq_zu),1.e-9), dq_zu ) 
     192      dt_zu = t_zu - T_s ;   dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 
     193      dq_zu = q_zu - q_s ;   dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 
    141194 
    142195      znu_a = visc_air(t_zt) ! Air viscosity (m^2/s) at zt given from temperature in (K) 
    143196 
    144       ztmp2 = 0.5 * 0.5  ! initial guess for wind gustiness contribution 
     197      ztmp2 = 0.5_wp*0.5_wp  ! initial guess for wind gustiness contribution 
    145198      U_blk = SQRT(U_zu*U_zu + ztmp2) 
    146199 
    147       ! z0     = 0.0001 
    148       ztmp2   = 10000.     ! optimization: ztmp2 == 1/z0 
     200      ztmp2   = 10000._wp     ! optimization: ztmp2 == 1/z0 (with z0 first guess == 0.0001) 
    149201      ztmp0   = LOG(zu*ztmp2) 
    150202      ztmp1   = LOG(10.*ztmp2) 
    151       u_star = 0.035*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10) 
    152  
    153       z0     = charn0*u_star*u_star/grav + 0.11*znu_a/u_star 
    154       z0t    = 0.1*EXP(vkarmn/(0.00115/(vkarmn/ztmp1)))   !  WARNING: 1/z0t ! 
    155  
    156       Cd     = (vkarmn/ztmp0)**2    ! first guess of Cd 
    157  
    158       ztmp0 = vkarmn*vkarmn/LOG(zt*z0t)/Cd 
     203      u_star = 0.035_wp*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10) 
     204 
     205      z0     = charn0*u_star*u_star/grav + 0.11_wp*znu_a/u_star 
     206      z0     = MIN(ABS(z0), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO 
     207      z0t    = 1._wp / ( 0.1_wp*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) ) 
     208      z0t    = MIN(ABS(z0t), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO 
     209 
     210      ztmp2  = vkarmn/ztmp0 
     211      Cd     = ztmp2*ztmp2    ! first guess of Cd 
     212 
     213      ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd 
    159214 
    160215      ztmp2 = Ri_bulk( zu, t_zu, dt_zu, q_zu, dq_zu, U_blk )   ! Ribu = Bulk Richardson number 
    161216 
    162217      !! First estimate of zeta_u, depending on the stability, ie sign of Ribu (ztmp2): 
    163       ztmp1 = 0.5 + SIGN( 0.5 , ztmp2 ) 
     218      ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 ) 
    164219      func_m = ztmp0*ztmp2 ! temporary array !! 
    165       !!             Ribu < 0                                 Ribu > 0   Beta = 1.25 
    166       func_h = (1.-ztmp1)*(func_m/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) &  ! temporary array !!! func_h == zeta_u 
    167          &  +     ztmp1*(func_m*(1. + 27./9.*ztmp2/ztmp0)) 
     220      func_h = (1.-ztmp1) * (func_m/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & !  Ribu < 0 ! temporary array !!! func_h == zeta_u 
     221         &  +     ztmp1   * (func_m*(1._wp + 27._wp/9._wp*ztmp2/func_m))              !  Ribu > 0 
     222      !#LB: should make sure that the "func_m" of "27./9.*ztmp2/func_m" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" ! 
    168223 
    169224      !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L 
    170       ztmp0   =        vkarmn/(LOG(zu*z0t) - psi_h_ecmwf(func_h)) 
     225      ztmp0  = vkarmn/(LOG(zu/z0t) - psi_h_ecmwf(func_h)) 
    171226 
    172227      u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_ecmwf(func_h)) 
     
    176231      ! What's need to be done if zt /= zu: 
    177232      IF( .NOT. l_zt_equal_zu ) THEN 
    178          ! 
    179233         !! First update of values at zu (or zt for wind) 
    180234         ztmp0 = psi_h_ecmwf(func_h) - psi_h_ecmwf(zt*func_h/zu)    ! zt*func_h/zu == zeta_t 
    181          ztmp1 = log(zt/zu) + ztmp0 
     235         ztmp1 = LOG(zt/zu) + ztmp0 
    182236         t_zu = t_zt - t_star/vkarmn*ztmp1 
    183237         q_zu = q_zt - q_star/vkarmn*ztmp1 
    184          q_zu = (0.5 + sign(0.5,q_zu))*q_zu !Makes it impossible to have negative humidity : 
    185  
    186          dt_zu = t_zu - sst  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu ) 
    187          dq_zu = q_zu - ssq  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu ) 
     238         q_zu = (0.5_wp + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity : 
    188239         ! 
    189       ENDIF 
     240         dt_zu = t_zu - T_s  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 
     241         dq_zu = q_zu - q_s  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 
     242      END IF 
    190243 
    191244 
     
    195248      !! First guess of inverse of Monin-Obukov length (1/L) : 
    196249      ztmp0 = (1. + rctv0*q_zu)  ! the factor to apply to temp. to get virt. temp... 
    197       Linv  =  grav*vkarmn*(t_star*ztmp0 + rctv0*t_zu*q_star) / ( u_star*u_star * t_zu*ztmp0 ) 
     250      Linv  =  grav*vkarmn*(t_star*ztmp0 + rctv0*t_zu*q_star) / MAX( u_star*u_star * t_zu*ztmp0 , 1.E-9 ) ! #LOLO 
     251      Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...) !#LOLO 
    198252 
    199253      !! Functions such as  u* = U_blk*vkarmn/func_m 
    200       ztmp1 = zu + z0 
    201       ztmp0 = ztmp1*Linv 
    202       func_m = LOG(ztmp1) -LOG(z0) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0*Linv) 
    203       func_h = LOG(ztmp1*z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(1./z0t*Linv) 
    204  
     254      ztmp0 = zu*Linv 
     255      func_m = LOG(zu) - LOG(z0)  - psi_m_ecmwf(ztmp0) + psi_m_ecmwf( z0*Linv) 
     256      func_h = LOG(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv) 
    205257 
    206258      !! ITERATION BLOCK 
    207       !! *************** 
    208  
    209259      DO j_itt = 1, nb_itt 
    210260 
     
    213263 
    214264         !! New estimate of the inverse of the Monin-Obukhon length (Linv == zeta/zu) : 
    215          Linv = ztmp0*func_m*func_m/func_h / zu     ! From Eq. 3.23, Chap.3, p.33, IFS doc - Cy31r1 
     265         Linv = ztmp0*func_m*func_m/func_h / zu     ! From Eq. 3.23, Chap.3.2.3, IFS doc - Cy40r1 
     266         !! Note: it is slightly different that the L we would get with the usual 
     267         !! expression, as in coare algorithm or in 'mod_phymbl.f90' (One_on_L_MO()) 
     268         Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...) !#LOLO 
    216269 
    217270         !! Update func_m with new Linv: 
    218          ztmp1 = zu + z0 
    219          func_m = LOG(ztmp1) -LOG(z0) - psi_m_ecmwf(ztmp1*Linv) + psi_m_ecmwf(z0*Linv) 
     271         func_m = LOG(zu) -LOG(z0) - psi_m_ecmwf(zu*Linv) + psi_m_ecmwf(z0*Linv) ! LB: should be "zu+z0" rather than "zu" alone, but z0 is tiny wrt zu! 
    220272 
    221273         !! Need to update roughness lengthes: 
     
    223275         ztmp2  = u_star*u_star 
    224276         ztmp1  = znu_a/u_star 
    225          z0    = alpha_M*ztmp1 + charn0*ztmp2/grav 
    226          z0t    = alpha_H*ztmp1                              ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1 
    227          z0q    = alpha_Q*ztmp1 
    228  
     277         z0     = MIN( ABS( alpha_M*ztmp1 + charn0*ztmp2/grav ) , 0.001_wp) 
     278         z0t    = MIN( ABS( alpha_H*ztmp1                     ) , 0.001_wp)   ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1 
     279         z0q    = MIN( ABS( alpha_Q*ztmp1                     ) , 0.001_wp) 
     280          
    229281         !! Update wind at 10m taking into acount convection-related wind gustiness: 
     282         !! => Chap. 3.2, IFS doc - Cy40r1, Eq.3.17 and Eq.3.18 + Eq.3.8 
    230283         ! Only true when unstable (L<0) => when ztmp0 < 0 => - !!! 
    231          ztmp2 = ztmp2 * (MAX(-zi0*Linv/vkarmn,0.))**(2./3.) ! => w*^2  (combining Eq. 3.8 and 3.18, hap.3, IFS doc - Cy31r1) 
     284         ztmp2 = ztmp2 * ( MAX(-zi0*Linv/vkarmn , 0._wp))**(2._wp/3._wp) ! => w*^2  (combining Eq. 3.8 and 3.18, hap.3, IFS doc - Cy31r1) 
    232285         !! => equivalent using Beta=1 (gustiness parameter, 1.25 for COARE, also zi0=600 in COARE..) 
    233          U_blk = MAX(sqrt(U_zu*U_zu + ztmp2), 0.2)              ! eq.3.17, Chap.3, p.32, IFS doc - Cy31r1 
     286         U_blk = MAX( SQRT(U_zu*U_zu + ztmp2) , 0.2_wp )              ! eq.3.17, Chap.3, p.32, IFS doc - Cy31r1 
    234287         ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0. 
    235288 
     
    238291         !! as well the air-sea differences: 
    239292         IF( .NOT. l_zt_equal_zu ) THEN 
    240  
    241293            !! Arrays func_m and func_h are free for a while so using them as temporary arrays... 
    242             func_h = psi_h_ecmwf((zu+z0)*Linv) ! temporary array !!! 
    243             func_m = psi_h_ecmwf((zt+z0)*Linv) ! temporary array !!! 
     294            func_h = psi_h_ecmwf(zu*Linv) ! temporary array !!! 
     295            func_m = psi_h_ecmwf(zt*Linv) ! temporary array !!! 
    244296 
    245297            ztmp2  = psi_h_ecmwf(z0t*Linv) 
    246298            ztmp0  = func_h - ztmp2 
    247             ztmp1  = vkarmn/(LOG(zu+z0) - LOG(z0t) - ztmp0) 
     299            ztmp1  = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0) 
    248300            t_star = dt_zu*ztmp1 
    249301            ztmp2  = ztmp0 - func_m + ztmp2 
     
    253305            ztmp2  = psi_h_ecmwf(z0q*Linv) 
    254306            ztmp0  = func_h - ztmp2 
    255             ztmp1  = vkarmn/(LOG(zu+z0) - LOG(z0q) - ztmp0) 
     307            ztmp1  = vkarmn/(LOG(zu) - LOG(z0q) - ztmp0) 
    256308            q_star = dq_zu*ztmp1 
    257309            ztmp2  = ztmp0 - func_m + ztmp2 
    258             ztmp1  = log(zt/zu) + ztmp2 
     310            ztmp1  = LOG(zt/zu) + ztmp2 
    259311            q_zu   = q_zt - q_star/vkarmn*ztmp1 
    260  
    261             dt_zu = t_zu - sst ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu ) 
    262             dq_zu = q_zu - ssq ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu ) 
    263  
    264312         END IF 
    265313 
    266314         !! Updating because of updated z0 and z0t and new Linv... 
    267          ztmp1 = zu + z0 
    268          ztmp0 = ztmp1*Linv 
    269          func_m = log(ztmp1) - LOG(z0 ) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0 *Linv) 
    270          func_h = log(ztmp1) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv) 
    271  
    272       END DO 
     315         ztmp0 = zu*Linv 
     316         func_m = log(zu) - LOG(z0 ) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0 *Linv) 
     317         func_h = log(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv) 
     318 
     319         !! SKIN related part 
     320         !! ----------------- 
     321         IF( l_use_skin ) THEN 
     322            !! compute transfer coefficients at zu : lolo: verifier... 
     323            Cd = vkarmn*vkarmn/(func_m*func_m) 
     324            Ch = vkarmn*vkarmn/(func_m*func_h) 
     325            ztmp1 = LOG(zu) - LOG(z0q) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0q*Linv)   ! func_q 
     326            Ce = vkarmn*vkarmn/(func_m*ztmp1) 
     327            ! Non-Solar heat flux to the ocean: 
     328            ztmp1 = U_blk*MAX(rho_air(t_zu, q_zu, slp), 1._wp)     ! rho*U10 
     329            ztmp2 = T_s*T_s 
     330            ztmp1 = ztmp1 * ( Ce*rLevap*(q_zu - q_s) + Ch*rCp_dry*(t_zu - T_s) ) & ! Total turb. heat flux 
     331               &       +    (rad_lw - emiss_w*stefan*ztmp2*ztmp2)                  ! Net longwave flux 
     332            !! Updating the values of the skin temperature T_s and q_s : 
     333            CALL CSWL_ECMWF( Qsw, ztmp1, u_star, zsst, T_s ) 
     334            q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp)  ! 200 -> just to avoid numerics problem on masked regions if silly values are given 
     335         END IF 
     336 
     337         IF( (l_use_skin).OR.(.NOT. l_zt_equal_zu) ) THEN 
     338            dt_zu = t_zu - T_s ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 
     339            dq_zu = q_zu - q_s ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 
     340         END IF 
     341 
     342      END DO !DO j_itt = 1, nb_itt 
    273343 
    274344      Cd = vkarmn*vkarmn/(func_m*func_m) 
    275345      Ch = vkarmn*vkarmn/(func_m*func_h) 
    276       ztmp1 = log((zu + z0)/z0q) - psi_h_ecmwf((zu + z0)*Linv) + psi_h_ecmwf(z0q*Linv)   ! func_q 
    277       Ce = vkarmn*vkarmn/(func_m*ztmp1) 
    278  
    279       ztmp1 = zu + z0 
    280       Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 )) 
    281       Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t)) 
    282       Cen = vkarmn*vkarmn / (log(ztmp1/z0q)*log(ztmp1/z0q)) 
     346      ztmp2 = log(zu/z0q) - psi_h_ecmwf(zu*Linv) + psi_h_ecmwf(z0q*Linv)   ! func_q 
     347      Ce = vkarmn*vkarmn/(func_m*ztmp2) 
     348 
     349      Cdn = vkarmn*vkarmn / (log(zu/z0 )*log(zu/z0 )) 
     350      Chn = vkarmn*vkarmn / (log(zu/z0t)*log(zu/z0t)) 
     351      Cen = vkarmn*vkarmn / (log(zu/z0q)*log(zu/z0q)) 
    283352 
    284353   END SUBROUTINE TURB_ECMWF 
     
    294363      !!         and L is M-O length 
    295364      !! 
    296       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     365      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    297366      !!---------------------------------------------------------------------------------- 
    298367      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ecmwf 
     
    306375         DO ji = 1, jpi 
    307376            ! 
    308             zzeta = MIN( pzeta(ji,jj) , 5. ) !! Very stable conditions (L positif and big!): 
     377            zzeta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!): 
    309378            ! 
    310379            ! Unstable (Paulson 1970): 
    311380            !   eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 
    312             zx = SQRT(ABS(1. - 16.*zzeta)) 
    313             ztmp = 1. + SQRT(zx) 
     381            zx = SQRT(ABS(1._wp - 16._wp*zzeta)) 
     382            ztmp = 1._wp + SQRT(zx) 
    314383            ztmp = ztmp*ztmp 
    315             psi_unst = LOG( 0.125*ztmp*(1. + zx) )   & 
    316                &       -2.*ATAN( SQRT(zx) ) + 0.5*rpi 
     384            psi_unst = LOG( 0.125_wp*ztmp*(1._wp + zx) )   & 
     385               &       -2._wp*ATAN( SQRT(zx) ) + 0.5_wp*rpi 
    317386            ! 
    318387            ! Unstable: 
    319388            ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 
    320             psi_stab = -2./3.*(zzeta - 5./0.35)*EXP(-0.35*zzeta) & 
    321                &       - zzeta - 2./3.*5./0.35 
     389            psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & 
     390               &       - zzeta - 2._wp/3._wp*5._wp/0.35_wp 
    322391            ! 
    323392            ! Combining: 
    324             stab = 0.5 + SIGN(0.5, zzeta) ! zzeta > 0 => stab = 1 
    325             ! 
    326             psi_m_ecmwf(ji,jj) = (1. - stab) * psi_unst & ! (zzeta < 0) Unstable 
    327                &                +      stab  * psi_stab   ! (zzeta > 0) Stable 
     393            stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 
     394            ! 
     395            psi_m_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable 
     396               &                +      stab  * psi_stab      ! (zzeta > 0) Stable 
    328397            ! 
    329398         END DO 
     
    332401   END FUNCTION psi_m_ecmwf 
    333402 
    334     
     403 
    335404   FUNCTION psi_h_ecmwf( pzeta ) 
    336405      !!---------------------------------------------------------------------------------- 
     
    342411      !!         and L is M-O length 
    343412      !! 
    344       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     413      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    345414      !!---------------------------------------------------------------------------------- 
    346415      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ecmwf 
     
    354423         DO ji = 1, jpi 
    355424            ! 
    356             zzeta = MIN(pzeta(ji,jj) , 5.)   ! Very stable conditions (L positif and big!): 
    357             ! 
    358             zx  = ABS(1. - 16.*zzeta)**.25        ! this is actually (1/phi_m)**2  !!! 
     425            zzeta = MIN(pzeta(ji,jj) , 5._wp)   ! Very stable conditions (L positif and big!): 
     426            ! 
     427            zx  = ABS(1._wp - 16._wp*zzeta)**.25        ! this is actually (1/phi_m)**2  !!! 
    359428            !                                     ! eq.3.19, Chap.3, p.33, IFS doc - Cy31r1 
    360429            ! Unstable (Paulson 1970) : 
    361             psi_unst = 2.*LOG(0.5*(1. + zx*zx))   ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 
     430            psi_unst = 2._wp*LOG(0.5_wp*(1._wp + zx*zx))   ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 
    362431            ! 
    363432            ! Stable: 
    364             psi_stab = -2./3.*(zzeta - 5./0.35)*EXP(-0.35*zzeta) & ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 
    365                &       - ABS(1. + 2./3.*zzeta)**1.5 - 2./3.*5./0.35 + 1.  
     433            psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 
     434               &       - ABS(1._wp + 2._wp/3._wp*zzeta)**1.5_wp - 2._wp/3._wp*5._wp/0.35_wp + 1._wp 
    366435            ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution... 
    367436            ! 
    368             stab = 0.5 + SIGN(0.5, zzeta) ! zzeta > 0 => stab = 1 
    369             ! 
    370             ! 
    371             psi_h_ecmwf(ji,jj) = (1. - stab) * psi_unst &   ! (zzeta < 0) Unstable 
    372                &                +    stab    * psi_stab     ! (zzeta > 0) Stable 
     437            stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 
     438            ! 
     439            ! 
     440            psi_h_ecmwf(ji,jj) = (1._wp - stab) * psi_unst &   ! (zzeta < 0) Unstable 
     441               &                +    stab    * psi_stab        ! (zzeta > 0) Stable 
    373442            ! 
    374443         END DO 
     
    382451      !! Bulk Richardson number (Eq. 3.25 IFS doc) 
    383452      !! 
    384       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     453      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    385454      !!---------------------------------------------------------------------------------- 
    386455      REAL(wp), DIMENSION(jpi,jpj) ::   Ri_bulk   ! 
     
    394463      !!---------------------------------------------------------------------------------- 
    395464      ! 
    396       Ri_bulk =   grav*pz/(pub*pub)                                          & 
    397          &      * ( pdt/(ptz - 0.5_wp*(pdt + grav*pz/(Cp_dry+Cp_vap*pqz)))  & 
     465      Ri_bulk =   grav*pz/(pub*pub)   & 
     466         &      * ( pdt/(ptz - 0.5_wp*(pdt + grav*pz/(rCp_dry + rCp_vap*pqz))) & 
    398467         &          + rctv0*pdq ) 
    399468      ! 
     
    401470 
    402471 
    403    FUNCTION visc_air(ptak) 
    404       !!---------------------------------------------------------------------------------- 
    405       !! Air kinetic viscosity (m^2/s) given from temperature in degrees... 
    406       !! 
    407       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    408       !!---------------------------------------------------------------------------------- 
    409       REAL(wp), DIMENSION(jpi,jpj)             ::   visc_air   ! 
    410       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak       ! air temperature in (K) 
    411       ! 
    412       INTEGER  ::   ji, jj      ! dummy loop indices 
    413       REAL(wp) ::   ztc, ztc2   ! local scalar 
    414       !!---------------------------------------------------------------------------------- 
    415       ! 
    416       DO jj = 1, jpj 
    417          DO ji = 1, jpi 
    418             ztc  = ptak(ji,jj) - rt0   ! air temp, in deg. C 
    419             ztc2 = ztc*ztc 
    420             visc_air(ji,jj) = 1.326e-5*(1. + 6.542E-3*ztc + 8.301e-6*ztc2 - 4.84e-9*ztc2*ztc) 
    421          END DO 
    422       END DO 
    423       ! 
    424    END FUNCTION visc_air 
    425472 
    426473   !!====================================================================== 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ncar.F90

    r10190 r11111  
    1111   !! 
    1212   !!       Routine turb_ncar maintained and developed in AeroBulk 
    13    !!                     (http://aerobulk.sourceforge.net/) 
     13   !!                     (https://github.com/brodeau/aerobulk/) 
    1414   !! 
    1515   !!                         L. Brodeau, 2015 
     
    4343 
    4444   PUBLIC ::   TURB_NCAR   ! called by sbcblk.F90 
    45  
    46    !                              ! NCAR own values for given constants: 
    47    REAL(wp), PARAMETER ::   rctv0 = 0.608   ! constant to obtain virtual temperature... 
    4845    
    4946   !!---------------------------------------------------------------------- 
     
    7673      !!      => 'vkarmn' and 'grav' 
    7774      !! 
    78       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     75      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    7976      !! 
    8077      !! INPUT : 
     
    238235      !! Origin: Large & Yeager 2008 eq.(11a) and eq.(11b) 
    239236      !! 
    240       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     237      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    241238      !!---------------------------------------------------------------------------------- 
    242239      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pw10           ! scalar wind speed at 10m (m/s) 
     
    277274      !!         and L is M-O length 
    278275      !! 
    279       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     276      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    280277      !!---------------------------------------------------------------------------------- 
    281278      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pzeta 
     
    311308      !!         and L is M-O length 
    312309      !! 
    313       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     310      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    314311      !!---------------------------------------------------------------------------------- 
    315312      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta 
Note: See TracChangeset for help on using the changeset viewer.