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 14592 – NEMO

Changeset 14592


Ignore:
Timestamp:
2021-03-05T17:03:57+01:00 (3 years ago)
Author:
gsamson
Message:

sbc_phy:

  • 2 new functions:
    • pres_temp: compute air pressure from air temperature and surface pressure
    • theta_exner: compute air/surface potential temperature following a dry adiabatic transformation
  • 3 modified functions:
    • Ri_bulk: use potential SST instead of SST + some simplifications TBD
    • bulk_formula: pass "rhoa" variable as input argument instead of calculing it with the old "gamma_moist" function inside the routine
    • update_qnsol_tau: pass "rhoa" variable as input argument to pass it to "bulk_formula" function
  • global small cleaning/cosmetics

sbcblk_algo_{ecmwf,coare3p?}: "rhoa" variable updated for CS/WL options and passed as input argument to "update_qnsol_tau" function

sbcblk & sbcabl:
add calls to "pres_temp" & "theta_exner" functions to compute air & surface potential temperatures (over both ocean and sea-ice) before computing sensible heat flux
special attention needed to check if CS/WL increment is properly applied

ticket #2632

Location:
NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/ABL/ablmod.F90

    r14433 r14592  
    104104      ! 
    105105      REAL(wp), DIMENSION(1:jpi,1:jpj       )    ::   zwnd_i, zwnd_j 
     106      REAL(wp), DIMENSION(1:jpi,1:jpj       )    ::   zsspt 
     107      REAL(wp), DIMENSION(1:jpi,1:jpj       )    ::   ztabs, zpre 
    106108      REAL(wp), DIMENSION(1:jpi      ,2:jpka)    ::   zCF 
    107109      ! 
    … …  
    142144      zrough(:,:) = z0_from_Cd( ght_abl(2), pCd_du(:,:) / MAX( pwndm(:,:), 0.5_wp ) ) ! #LB: z0_from_Cd is define in sbc_phy.F90... 
    143145 
     146      ! sea surface potential temperature [K] 
     147      zsspt(:,:) = theta_exner( psst(:,:)+rt0, pslp_dta(:,:) ) 
     148 
    144149      ! 
    145150      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    … …  
    187192               DO ji = 1,jpi  ! surface boundary condition for temperature 
    188193                  zztmp1 = psen(ji, jj) 
    189                   zztmp2 = psen(ji, jj) * ( psst(ji, jj) + rt0 ) 
     194                  zztmp2 = psen(ji, jj) * zsspt(ji, jj) 
    190195#if defined key_si3 
    191196                  zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj) 
    … …  
    495500               zcff2 = jp_alp3_dyn * zsig**3 + jp_alp2_dyn * zsig**2   & 
    496501                  &  + jp_alp1_dyn * zsig    + jp_alp0_dyn 
    497                zcff  = (1._wp-zmsk) + zmsk * zcff2 * rn_Dt   ! zcff = 1 for masked points 
    498                                                              ! rn_Dt = rDt_abl / nn_fsbc 
     502               zcff  = (1._wp-zmsk) + zmsk * zcff2 * rDt_abl   ! zcff = 1 for masked points 
    499503               zcff  = zcff * rest_eq(ji,jj) 
    500504               u_abl( ji, jj, jk, nt_a ) = (1._wp - zcff ) *  u_abl( ji, jj, jk, nt_a )   & 
    … …  
    518522            zcff2 = jp_alp3_tra * zsig**3 + jp_alp2_tra * zsig**2   & 
    519523               &  + jp_alp1_tra * zsig    + jp_alp0_tra 
    520             zcff  = (1._wp-zmsk) + zmsk * zcff2 * rn_Dt   ! zcff = 1 for masked points 
    521                                                           ! rn_Dt = rDt_abl / nn_fsbc 
     524            zcff  = (1._wp-zmsk) + zmsk * zcff2 * rDt_abl   ! zcff = 1 for masked points 
    522525            tq_abl( ji, jj, jk, nt_a, jp_ta ) = (1._wp - zcff ) * tq_abl( ji, jj, jk, nt_a, jp_ta )   & 
    523526               &                                       + zcff   * pt_dta( ji, jj, jk ) 
    … …  
    586589      ! 
    587590      DO_2D( 1, 1, 1, 1 ) 
    588          ztemp          =  tq_abl( ji, jj, 2, nt_a, jp_ta ) 
    589          zhumi          =  tq_abl( ji, jj, 2, nt_a, jp_qa ) 
    590          zcff           = rho_air( ztemp, zhumi, pslp_dta( ji, jj ) ) 
    591          psen( ji, jj ) =    - cp_air(zhumi) * zcff * psen(ji,jj) * ( psst(ji,jj) + rt0 - ztemp )   !GS: negative sign to respect aerobulk convention 
    592          pevp( ji, jj ) = rn_efac*MAX( 0._wp,  zcff * pevp(ji,jj) * ( pssq(ji,jj)       - zhumi ) ) 
     591         ztemp          = tq_abl( ji, jj, 2, nt_a, jp_ta ) 
     592         zhumi          = tq_abl( ji, jj, 2, nt_a, jp_qa ) 
     593         zpre( ji, jj ) = pres_temp( zhumi, pslp_dta(ji,jj), ght_abl(2), ptpot=ztemp, pta=ztabs( ji, jj ) ) 
     594         zcff           = rho_air( ztabs( ji, jj ), zhumi, zpre( ji, jj ) ) 
     595         psen( ji, jj ) =    - cp_air(zhumi) * zcff * psen(ji,jj) * ( zsspt(ji,jj) - ztemp )         !GS: negative sign to respect aerobulk convention 
     596         pevp( ji, jj ) = rn_efac*MAX( 0._wp,  zcff * pevp(ji,jj) * ( pssq(ji,jj)  - zhumi ) ) 
    593597         plat( ji, jj ) = - L_vap( psst(ji,jj) ) * pevp( ji, jj ) 
    594598         rhoa( ji, jj ) = zcff 
    595599      END_2D 
     600      !!GS: debug only; to be removed 
     601      CALL iom_put ( "pres_zu", zpre(:,:) )  
     602      CALL iom_put ( "tabs_zu", ztabs(:,:) )  
    596603 
    597604      DO_2D( 0, 1, 0, 1 ) 
    … …  
    639646      !    Wind stress relative to the moving ice ( U10m - U_ice )   ! 
    640647      ! ------------------------------------------------------------ ! 
    641       DO_2D( 0, 0, 0, 0 ) 
    642          ptaui_ice(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) + rhoa(ji,jj) * pCd_du_ice(ji,jj)      )   & 
    643             &                      * ( 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) ) - pssu_ice(ji,jj) ) 
    644          ptauj_ice(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) + rhoa(ji,jj) * pCd_du_ice(ji,jj)      )   & 
    645             &                      * ( 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) - pssv_ice(ji,jj) ) 
    646       END_2D 
    647       CALL lbc_lnk( 'ablmod', ptaui_ice, 'U', -1.0_wp, ptauj_ice, 'V', -1.0_wp ) 
    648       ! 
    649       IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=ptaui_ice  , clinfo1=' abl_stp: putaui : '   & 
    650          &                                , tab2d_2=ptauj_ice  , clinfo2='          pvtaui : ' ) 
     648      !DO_2D( 0, 0, 0, 0 ) 
     649      !   ptaui_ice(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) + rhoa(ji,jj) * pCd_du_ice(ji,jj)      )   & 
     650      !      &                      * ( 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) ) - pssu_ice(ji,jj) ) 
     651      !   ptauj_ice(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) + rhoa(ji,jj) * pCd_du_ice(ji,jj)      )   & 
     652      !      &                      * ( 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) - pssv_ice(ji,jj) ) 
     653      !END_2D 
     654      !CALL lbc_lnk( 'ablmod', ptaui_ice, 'U', -1.0_wp, ptauj_ice, 'V', -1.0_wp ) 
     655      !! 
     656      !IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=ptaui_ice  , clinfo1=' abl_stp: putaui : '   & 
     657      !   &                                , tab2d_2=ptauj_ice  , clinfo2='          pvtaui : ' ) 
     658 
    651659      ! ------------------------------------------------------------ ! 
    652660      !    Wind stress relative to the moving ice ( U10m - U_ice )   ! 

  • NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/ABL/sbcabl.F90

    r14239 r14592  
    218218      IF(lwp) THEN 
    219219         IF(nn_dyn_restore > 0) THEN 
    220             WRITE(numout,*) ' ABL Minimum value for dynamics restoring = ',zcff 
    221             WRITE(numout,*) ' ABL Maximum value for dynamics restoring = ',zcff1 
     220            WRITE(numout,*) ' Minimum value for ABL dynamics restoring = ',zcff 
     221            WRITE(numout,*) ' Maximum value for ABL dynamics restoring = ',zcff1 
    222222            ! Check that restoring coefficients are between 0 and 1 
    223223            IF( zcff1 - nn_fsbc > 0.001_wp .OR. zcff1 < 0._wp )   & 
    … …  
    236236         &    + jp_alp1_tra * jp_bmax    + jp_alp0_tra ) * rDt_abl 
    237237      IF(lwp) THEN 
    238          WRITE(numout,*) ' ABL Minimum value for tracers restoring = ',zcff 
    239          WRITE(numout,*) ' ABL Maximum value for tracers restoring = ',zcff1 
     238         WRITE(numout,*) ' Minimum value for ABL tracers restoring = ',zcff 
     239         WRITE(numout,*) ' Maximum value for ABL tracers restoring = ',zcff1 
    240240         ! Check that restoring coefficients are between 0 and 1 
    241241         IF( zcff1 - nn_fsbc > 0.001_wp .OR. zcff1 < 0._wp )   & 

  • NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/OCE/SBC/sbc_phy.F90

    r14110 r14592  
    2323 
    2424   IMPLICIT NONE 
    25    !PRIVATE 
    26    PUBLIC !! Haleluja that was the solution 
     25   PUBLIC !! Haleluja that was the solution for AGRIF 
    2726 
    2827   INTEGER , PARAMETER, PUBLIC  :: nb_iter0 = 5 ! Default number of itterations in bulk-param algorithms (can be overriden b.m.o `nb_iter` optional argument) 
    2928 
    3029   !!  (mainly removed from sbcblk.F90) 
    31    REAL(wp), PARAMETER, PUBLIC :: rCp_dry = 1005.0_wp   !: Specic heat of dry air, constant pressure      [J/K/kg] 
    32    REAL(wp), PARAMETER, PUBLIC :: rCp_vap = 1860.0_wp   !: Specic heat of water vapor, constant pressure  [J/K/kg] 
    33    REAL(wp), PARAMETER, PUBLIC :: R_dry   = 287.05_wp   !: Specific gas constant for dry air              [J/K/kg] 
    34    REAL(wp), PARAMETER, PUBLIC :: R_vap   = 461.495_wp  !: Specific gas constant for water vapor          [J/K/kg] 
    35    REAL(wp), PARAMETER, PUBLIC :: reps0   = R_dry/R_vap !: ratio of gas constant for dry air and water vapor => ~ 0.622 
    36    REAL(wp), PARAMETER, PUBLIC :: rctv0   = R_vap/R_dry - 1._wp  !: for virtual temperature (== (1-eps)/eps) => ~ 0.608 
    37    REAL(wp), PARAMETER, PUBLIC :: rCp_air = 1000.5_wp   !: specific heat of air (only used for ice fluxes now...) 
    38    REAL(wp), PARAMETER, PUBLIC :: albo    = 0.066_wp    !: ocean albedo assumed to be constant 
     30   REAL(wp), PARAMETER, PUBLIC :: rCp_dry = 1005.0_wp             !: Specic heat of dry air, constant pressure       [J/K/kg] 
     31   REAL(wp), PARAMETER, PUBLIC :: rCp_vap = 1860.0_wp             !: Specic heat of water vapor, constant pressure   [J/K/kg] 
     32   REAL(wp), PARAMETER, PUBLIC :: R_dry   = 287.05_wp             !: Specific gas constant for dry air               [J/K/kg] 
     33   REAL(wp), PARAMETER, PUBLIC :: R_vap   = 461.495_wp            !: Specific gas constant for water vapor           [J/K/kg] 
     34   REAL(wp), PARAMETER, PUBLIC :: reps0   = R_dry/R_vap           !: ratio of gas constant for dry air and water vapor => ~ 0.622 
     35   REAL(wp), PARAMETER, PUBLIC :: rctv0   = R_vap/R_dry - 1._wp   !: for virtual temperature (== (1-eps)/eps) => ~ 0.608 
     36   REAL(wp), PARAMETER, PUBLIC :: rCp_air = 1000.5_wp             !: specific heat of air (only used for ice fluxes now...) 
     37   REAL(wp), PARAMETER, PUBLIC :: albo    = 0.066_wp              !: ocean albedo assumed to be constant 
     38   REAL(wp), PARAMETER, PUBLIC :: R_gas   = 8.314510_wp           !: Universal molar gas constant                    [J/mol/K]  
     39   REAL(wp), PARAMETER, PUBLIC :: rmm_dryair = 28.9647e-3_wp      !: dry air molar mass / molecular weight           [kg/mol] 
     40   REAL(wp), PARAMETER, PUBLIC :: rmm_water  = 18.0153e-3_wp      !: water   molar mass / molecular weight           [kg/mol] 
     41   REAL(wp), PARAMETER, PUBLIC :: rmm_ratio  = rmm_water / rmm_dryair 
     42   REAL(wp), PARAMETER, PUBLIC :: rgamma_dry = R_gas / ( rmm_dryair * rCp_dry )  !: Poisson constant for dry air 
     43   REAL(wp), PARAMETER, PUBLIC :: rpref      = 1.e5_wp            !: reference air pressure for exner function       [Pa] 
    3944   ! 
    4045   REAL(wp), PARAMETER, PUBLIC :: rho0_a  = 1.2_wp      !: Approx. of density of air                       [kg/m^3] 
    … …  
    8388   END INTERFACE virt_temp 
    8489 
     90   INTERFACE pres_temp 
     91      MODULE PROCEDURE pres_temp_vctr, pres_temp_sclr 
     92   END INTERFACE pres_temp 
     93 
     94   INTERFACE theta_exner 
     95      MODULE PROCEDURE theta_exner_vctr, theta_exner_sclr 
     96   END INTERFACE theta_exner 
     97 
    8598   INTERFACE visc_air 
    8699      MODULE PROCEDURE visc_air_vctr, visc_air_sclr 
    … …  
    98111      MODULE PROCEDURE e_sat_ice_vctr, e_sat_ice_sclr 
    99112   END INTERFACE e_sat_ice 
     113 
    100114   INTERFACE de_sat_dt_ice 
    101115      MODULE PROCEDURE de_sat_dt_ice_vctr, de_sat_dt_ice_sclr 
    … …  
    148162 
    149163   PUBLIC virt_temp 
     164   PUBLIC pres_temp 
     165   PUBLIC theta_exner 
    150166   PUBLIC rho_air 
    151167   PUBLIC visc_air 
    … …  
    158174   PUBLIC q_air_rh 
    159175   PUBLIC dq_sat_dt_ice 
    160    !: 
     176   ! 
    161177   PUBLIC update_qnsol_tau 
    162178   PUBLIC alpha_sw 
    … …  
    205221      ! 
    206222   END FUNCTION virt_temp_sclr 
    207    !! 
     223 
    208224   FUNCTION virt_temp_vctr( pta, pqa ) 
     225 
    209226      REAL(wp), DIMENSION(jpi,jpj)             :: virt_temp_vctr !: virtual temperature [K] 
    210227      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta !: absolute or potential air temperature [K] 
    211228      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa !: specific humidity of air   [kg/kg] 
     229 
    212230      virt_temp_vctr(:,:) = pta(:,:) * (1._wp + rctv0*pqa(:,:)) 
     231 
    213232   END FUNCTION virt_temp_vctr 
    214    !=============================================================================================== 
     233 
     234 
     235   FUNCTION pres_temp_sclr( pqspe, pslp, pz, ptpot, pta, l_ice ) 
     236 
     237      !!------------------------------------------------------------------------------- 
     238      !!                           ***  FUNCTION pres_temp  *** 
     239      !! 
     240      !! ** Purpose : compute air pressure using barometric equation 
     241      !!              from either potential or absolute air temperature 
     242      !! ** Author: G. Samson, Feb 2021 
     243      !!------------------------------------------------------------------------------- 
     244 
     245      REAL(wp)                          :: pres_temp_sclr    ! air pressure              [Pa] 
     246      REAL(wp), INTENT(in )             :: pqspe             ! air specific humidity     [kg/kg] 
     247      REAL(wp), INTENT(in )             :: pslp              ! sea-level pressure        [Pa] 
     248      REAL(wp), INTENT(in )             :: pz                ! height above surface      [m] 
     249      REAL(wp), INTENT(in )  , OPTIONAL :: ptpot             ! air potential temperature [K] 
     250      REAL(wp), INTENT(inout), OPTIONAL :: pta               ! air absolute temperature  [K] 
     251      REAL(wp)                          :: ztpot, zta, zpa, zxm, zmask, zqsat 
     252      INTEGER                           :: it, niter = 3     ! iteration indice and number 
     253      LOGICAL , INTENT(in)   , OPTIONAL :: l_ice             ! sea-ice presence 
     254      LOGICAL                           :: lice              ! sea-ice presence 
     255 
     256      IF( PRESENT(ptpot) ) THEN 
     257        zmask = 1._wp 
     258        ztpot = ptpot 
     259      ELSE 
     260        zmask = 0._wp  
     261        zta   = pta 
     262      ENDIF 
     263 
     264      lice = .FALSE. 
     265      IF( PRESENT(l_ice) ) lice = l_ice 
     266 
     267      zpa = pslp              ! air pressure first guess [Pa] 
     268      DO it = 1, niter 
     269         zta   = ztpot * ( zpa / rpref )**rgamma_dry * zmask + (1._wp - zmask) * zta 
     270         zqsat = q_sat( zta, zpa, l_ice=lice )                                   ! saturation specific humidity [kg/kg] 
     271         zxm   = (1._wp - pqspe/zqsat) * rmm_dryair + pqspe/zqsat * rmm_water    ! moist air molar mass [kg/mol] 
     272         zpa   = pslp * EXP( -grav * zxm * pz / ( R_gas * zta ) ) 
     273      END DO 
     274 
     275      pres_temp_sclr = zpa 
     276      IF(( PRESENT(pta) ).AND.( PRESENT(ptpot) )) pta = zta 
     277 
     278   END FUNCTION pres_temp_sclr 
     279 
     280 
     281   FUNCTION pres_temp_vctr( pqspe, pslp, pz, ptpot, pta, l_ice ) 
     282 
     283      !!------------------------------------------------------------------------------- 
     284      !!                           ***  FUNCTION pres_temp  *** 
     285      !! 
     286      !! ** Purpose : compute air pressure using barometric equation 
     287      !!              from either potential or absolute air temperature 
     288      !! ** Author: G. Samson, Feb 2021 
     289      !!------------------------------------------------------------------------------- 
     290 
     291      REAL(wp), DIMENSION(jpi,jpj)                          :: pres_temp_vctr    ! air pressure              [Pa] 
     292      REAL(wp), DIMENSION(jpi,jpj), INTENT(in )             :: pqspe             ! air specific humidity     [kg/kg] 
     293      REAL(wp), DIMENSION(jpi,jpj), INTENT(in )             :: pslp              ! sea-level pressure        [Pa] 
     294      REAL(wp),                     INTENT(in )             :: pz                ! height above surface      [m] 
     295      REAL(wp), DIMENSION(jpi,jpj), INTENT(in )  , OPTIONAL :: ptpot             ! air potential temperature [K] 
     296      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout), OPTIONAL :: pta               ! air absolute temperature  [K] 
     297      INTEGER                                               :: ji, jj            ! loop indices 
     298      LOGICAL                     , INTENT(in)   , OPTIONAL :: l_ice             ! sea-ice presence 
     299      LOGICAL                                               :: lice              ! sea-ice presence 
     300  
     301      lice = .FALSE. 
     302      IF( PRESENT(l_ice) ) lice = l_ice 
     303 
     304      IF( PRESENT(ptpot) ) THEN 
     305         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     306            pres_temp_vctr(ji,jj) = pres_temp_sclr( pqspe(ji,jj), pslp(ji,jj), pz, ptpot=ptpot(ji,jj), pta=pta(ji,jj), l_ice=lice ) 
     307         END_2D 
     308      ELSE 
     309         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     310            pres_temp_vctr(ji,jj) = pres_temp_sclr( pqspe(ji,jj), pslp(ji,jj), pz, pta=pta(ji,jj), l_ice=lice ) 
     311         END_2D 
     312      ENDIF 
     313 
     314   END FUNCTION pres_temp_vctr 
     315 
     316 
     317   FUNCTION theta_exner_sclr( pta, ppa ) 
     318 
     319      !!------------------------------------------------------------------------------- 
     320      !!                           ***  FUNCTION theta_exner  *** 
     321      !! 
     322      !! ** Purpose : compute air/surface potential temperature from absolute temperature 
     323      !!              and pressure using Exner function 
     324      !! ** Author: G. Samson, Feb 2021 
     325      !!------------------------------------------------------------------------------- 
     326 
     327      REAL(wp)             :: theta_exner_sclr   ! air/surface potential temperature [K] 
     328      REAL(wp), INTENT(in) :: pta                ! air/surface absolute temperature  [K] 
     329      REAL(wp), INTENT(in) :: ppa                ! air/surface pressure              [Pa] 
     330       
     331      theta_exner_sclr = pta * ( rpref / ppa ) ** rgamma_dry 
     332 
     333   END FUNCTION theta_exner_sclr 
     334 
     335   FUNCTION theta_exner_vctr( pta, ppa ) 
     336 
     337      !!------------------------------------------------------------------------------- 
     338      !!                           ***  FUNCTION theta_exner  *** 
     339      !! 
     340      !! ** Purpose : compute air/surface potential temperature from absolute temperature 
     341      !!              and pressure using Exner function 
     342      !! ** Author: G. Samson, Feb 2021 
     343      !!------------------------------------------------------------------------------- 
     344 
     345      REAL(wp), DIMENSION(jpi,jpj)             :: theta_exner_vctr   ! air/surface potential temperature [K] 
     346      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta                ! air/surface absolute temperature  [K] 
     347      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa                ! air/surface pressure              [Pa] 
     348      INTEGER                                  :: ji, jj             ! loop indices 
     349 
     350      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     351         theta_exner_vctr(ji,jj) = theta_exner_sclr( pta(ji,jj), ppa(ji,jj) ) 
     352      END_2D 
     353 
     354   END FUNCTION theta_exner_vctr 
    215355 
    216356 
    … …  
    228368      REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air_vctr   ! density of moist air   [kg/m^3] 
    229369      !!------------------------------------------------------------------------------- 
     370 
    230371      rho_air_vctr = MAX( ppa / (R_dry*ptak * ( 1._wp + rctv0*pqa )) , 0.8_wp ) 
     372 
    231373   END FUNCTION rho_air_vctr 
    232374 
    … …  
    245387      !!------------------------------------------------------------------------------- 
    246388      rho_air_sclr = MAX( ppa / (R_dry*ptak * ( 1._wp + rctv0*pqa )) , 0.8_wp ) 
     389 
    247390   END FUNCTION rho_air_sclr 
    248  
    249  
    250391 
    251392 
    … …  
    269410 
    270411   FUNCTION visc_air_vctr(ptak) 
     412 
    271413      REAL(wp), DIMENSION(jpi,jpj)             ::   visc_air_vctr   ! kinetic viscosity (m^2/s) 
    272414      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak       ! air temperature in (K) 
    273415      INTEGER  ::   ji, jj      ! dummy loop indices 
    274       DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    275       visc_air_vctr(ji,jj) = visc_air_sclr( ptak(ji,jj) ) 
    276       END_2D 
     416 
     417      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     418         visc_air_vctr(ji,jj) = visc_air_sclr( ptak(ji,jj) ) 
     419      END_2D 
     420 
    277421   END FUNCTION visc_air_vctr 
    278422 
    … …  
    319463      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    320464      !!------------------------------------------------------------------------------- 
    321       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg] 
     465      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa           ! air specific humidity         [kg/kg] 
    322466      REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air_vctr   ! specific heat of moist air   [J/K/kg] 
    323467      !!------------------------------------------------------------------------------- 
     468 
    324469      cp_air_vctr = rCp_dry + rCp_vap * pqa 
     470 
    325471   END FUNCTION cp_air_vctr 
    326472 
    … …  
    336482      REAL(wp)             :: cp_air_sclr   ! specific heat of moist air   [J/K/kg] 
    337483      !!------------------------------------------------------------------------------- 
     484 
    338485      cp_air_sclr = rCp_dry + rCp_vap * pqa 
     486 
    339487   END FUNCTION cp_air_sclr 
    340488 
    341489 
    342  
    343    !=============================================================================================== 
    344490   FUNCTION gamma_moist_sclr( ptak, pqa ) 
    345491      !!---------------------------------------------------------------------------------- 
    … …  
    366512      !! 
    367513   END FUNCTION gamma_moist_sclr 
    368    !! 
     514 
    369515   FUNCTION gamma_moist_vctr( ptak, pqa ) 
     516 
    370517      REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist_vctr 
    371518      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak 
    372519      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa 
    373520      INTEGER  :: ji, jj 
    374       DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    375       gamma_moist_vctr(ji,jj) = gamma_moist_sclr( ptak(ji,jj), pqa(ji,jj) ) 
    376       END_2D 
     521 
     522      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     523         gamma_moist_vctr(ji,jj) = gamma_moist_sclr( ptak(ji,jj), pqa(ji,jj) ) 
     524      END_2D 
     525 
    377526   END FUNCTION gamma_moist_vctr 
    378    !=============================================================================================== 
    379527 
    380528 
    … …  
    399547      ! 
    400548      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    401       ! 
    402       zqa = (1._wp + rctv0*pqa(ji,jj)) 
    403       ! 
    404       ! The main concern is to know whether, the vertical turbulent flux of virtual temperature, < u' theta_v' > is estimated with: 
    405       !  a/  -u* [ theta* (1 + 0.61 q) + 0.61 theta q* ] => this is the one that seems correct! chose this one! 
    406       !                      or 
    407       !  b/  -u* [ theta*              + 0.61 theta q* ] 
    408       ! 
    409       One_on_L(ji,jj) = grav*vkarmn*( pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj) ) & 
    410          &               / MAX( pus(ji,jj)*pus(ji,jj)*ptha(ji,jj)*zqa , 1.E-9_wp ) 
    411       ! 
     549         zqa = (1._wp + rctv0*pqa(ji,jj)) 
     550         ! 
     551         ! The main concern is to know whether, the vertical turbulent flux of virtual temperature, < u' theta_v' > is estimated with: 
     552         !  a/  -u* [ theta* (1 + 0.61 q) + 0.61 theta q* ] => this is the one that seems correct! chose this one! 
     553         !                      or 
     554         !  b/  -u* [ theta*              + 0.61 theta q* ] 
     555         ! 
     556         One_on_L(ji,jj) = grav*vkarmn*( pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj) ) & 
     557            &               / MAX( pus(ji,jj)*pus(ji,jj)*ptha(ji,jj)*zqa , 1.E-9_wp ) 
    412558      END_2D 
    413559      ! 
    … …  
    417563 
    418564 
    419    !=============================================================================================== 
    420565   FUNCTION Ri_bulk_sclr( pz, psst, ptha, pssq, pqa, pub,  pta_layer, pqa_layer ) 
    421566      !!---------------------------------------------------------------------------------- 
    … …  
    428573      REAL(wp)             :: Ri_bulk_sclr 
    429574      REAL(wp), INTENT(in) :: pz    ! height above the sea (aka "delta z")  [m] 
    430       REAL(wp), INTENT(in) :: psst  ! SST                                   [K] 
     575      REAL(wp), INTENT(in) :: psst  ! potential SST                         [K] 
    431576      REAL(wp), INTENT(in) :: ptha  ! pot. air temp. at height "pz"         [K] 
    432577      REAL(wp), INTENT(in) :: pssq  ! 0.98*q_sat(SST)                   [kg/kg] 
    … …  
    438583      LOGICAL  :: l_ptqa_l_prvd = .FALSE. 
    439584      REAL(wp) :: zqa, zta, zgamma, zdthv, ztv, zsstv  ! local scalars 
     585      REAL(wp) :: ztptv 
    440586      !!------------------------------------------------------------------- 
    441587      IF( PRESENT(pta_layer) .AND. PRESENT(pqa_layer) ) l_ptqa_l_prvd=.TRUE. 
    442588      ! 
    443       zsstv = virt_temp_sclr( psst, pssq )          ! virtual SST (absolute==potential because z=0!) 
    444       ! 
    445       zdthv = virt_temp_sclr( ptha, pqa  ) - zsstv  ! air-sea delta of "virtual potential temperature" 
    446       ! 
    447       !! ztv: estimate of the ABSOLUTE virtual temp. within the layer 
    448       IF( l_ptqa_l_prvd ) THEN 
    449          ztv = virt_temp_sclr( pta_layer, pqa_layer ) 
    450       ELSE 
    451          zqa = 0.5_wp*( pqa  + pssq )                             ! ~ mean q within the layer... 
    452          zta = 0.5_wp*( psst + ptha - gamma_moist(ptha, zqa)*pz ) ! ~ mean absolute temperature of air within the layer 
    453          zta = 0.5_wp*( psst + ptha - gamma_moist( zta, zqa)*pz ) ! ~ mean absolute temperature of air within the layer 
    454          zgamma =  gamma_moist(zta, zqa)                          ! Adiabatic lapse-rate for moist air within the layer 
    455          ztv = 0.5_wp*( zsstv + virt_temp_sclr( ptha-zgamma*pz, pqa ) ) 
    456       END IF 
    457       ! 
    458       Ri_bulk_sclr = grav*zdthv*pz / ( ztv*pub*pub )      ! the usual definition of Ri_bulk_sclr 
     589      zsstv = virt_temp_sclr( psst, pssq )   ! virtual potential SST 
     590      ztptv = virt_temp_sclr( ptha, pqa  )   ! virtual potential air temperature 
     591      zdthv = ztptv - zsstv                  ! air-sea delta of "virtual potential temperature" 
     592      ! 
     593      Ri_bulk_sclr = grav * zdthv * pz / ( ztptv * pub * pub )      ! the usual definition of Ri_bulk_sclr 
    459594      ! 
    460595   END FUNCTION Ri_bulk_sclr 
    461    !! 
     596 
    462597   FUNCTION Ri_bulk_vctr( pz, psst, ptha, pssq, pqa, pub,  pta_layer, pqa_layer ) 
     598 
    463599      REAL(wp), DIMENSION(jpi,jpj)             :: Ri_bulk_vctr 
    464600      REAL(wp)                    , INTENT(in) :: pz    ! height above the sea (aka "delta z")  [m] 
    … …  
    473609      LOGICAL  :: l_ptqa_l_prvd = .FALSE. 
    474610      INTEGER  ::   ji, jj 
     611 
    475612      IF( PRESENT(pta_layer) .AND. PRESENT(pqa_layer) ) l_ptqa_l_prvd=.TRUE. 
    476613      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    477       IF( l_ptqa_l_prvd ) THEN 
    478          Ri_bulk_vctr(ji,jj) = Ri_bulk_sclr( pz, psst(ji,jj), ptha(ji,jj), pssq(ji,jj), pqa(ji,jj), pub(ji,jj), & 
    479             &                                pta_layer=pta_layer(ji,jj ),  pqa_layer=pqa_layer(ji,jj ) ) 
    480       ELSE 
    481          Ri_bulk_vctr(ji,jj) = Ri_bulk_sclr( pz, psst(ji,jj), ptha(ji,jj), pssq(ji,jj), pqa(ji,jj), pub(ji,jj) ) 
    482       END IF 
    483       END_2D 
     614         IF( l_ptqa_l_prvd ) THEN  !!GS: "IF" inside loop needs to be removed 
     615            Ri_bulk_vctr(ji,jj) = Ri_bulk_sclr( pz, psst(ji,jj), ptha(ji,jj), pssq(ji,jj), pqa(ji,jj), pub(ji,jj), & 
     616               &                                pta_layer=pta_layer(ji,jj ),  pqa_layer=pqa_layer(ji,jj ) ) 
     617         ELSE 
     618            Ri_bulk_vctr(ji,jj) = Ri_bulk_sclr( pz, psst(ji,jj), ptha(ji,jj), pssq(ji,jj), pqa(ji,jj), pub(ji,jj) ) 
     619         END IF 
     620      END_2D 
     621 
    484622   END FUNCTION Ri_bulk_vctr 
    485    !=============================================================================================== 
    486  
    487    !=============================================================================================== 
     623 
     624 
    488625   FUNCTION e_sat_sclr( ptak ) 
    489626      !!---------------------------------------------------------------------------------- 
    … …  
    510647      ! 
    511648   END FUNCTION e_sat_sclr 
    512    !! 
     649 
    513650   FUNCTION e_sat_vctr(ptak) 
    514651      REAL(wp), DIMENSION(jpi,jpj)             :: e_sat_vctr !: vapour pressure at saturation  [Pa] 
    … …  
    519656      END_2D 
    520657   END FUNCTION e_sat_vctr 
    521    !=============================================================================================== 
    522  
    523    !=============================================================================================== 
     658 
     659 
    524660   FUNCTION e_sat_ice_sclr(ptak) 
    525661      !!--------------------------------------------------------------------------------- 
    … …  
    537673      !! 
    538674      e_sat_ice_sclr = 100._wp * 10._wp**zle 
     675    
    539676   END FUNCTION e_sat_ice_sclr 
    540    !! 
     677 
    541678   FUNCTION e_sat_ice_vctr(ptak) 
    542679      !! Same as "e_sat" but over ice rather than water! 
    … …  
    546683      !!---------------------------------------------------------------------------------- 
    547684      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    548       e_sat_ice_vctr(ji,jj) = e_sat_ice_sclr( ptak(ji,jj) ) 
    549       END_2D 
     685         e_sat_ice_vctr(ji,jj) = e_sat_ice_sclr( ptak(ji,jj) ) 
     686      END_2D 
     687 
    550688   END FUNCTION e_sat_ice_vctr 
    551    !! 
     689 
     690 
    552691   FUNCTION de_sat_dt_ice_sclr(ptak) 
    553692      !!--------------------------------------------------------------------------------- 
    … …  
    567706      de_sat_dt_ice_sclr = LOG(10._wp) * zde * e_sat_ice_sclr(zta) 
    568707   END FUNCTION de_sat_dt_ice_sclr 
    569    !! 
     708 
    570709   FUNCTION de_sat_dt_ice_vctr(ptak) 
    571710      !! Same as "e_sat" but over ice rather than water! 
    … …  
    575714      !!---------------------------------------------------------------------------------- 
    576715      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    577       de_sat_dt_ice_vctr(ji,jj) = de_sat_dt_ice_sclr( ptak(ji,jj) ) 
    578       END_2D 
     716         de_sat_dt_ice_vctr(ji,jj) = de_sat_dt_ice_sclr( ptak(ji,jj) ) 
     717      END_2D 
     718 
    579719   END FUNCTION de_sat_dt_ice_vctr 
    580720 
    581721 
    582  
    583    !=============================================================================================== 
    584  
    585    !=============================================================================================== 
    586722   FUNCTION q_sat_sclr( pta, ppa,  l_ice ) 
    587723      !!--------------------------------------------------------------------------------- 
    … …  
    607743      END IF 
    608744      q_sat_sclr = reps0*ze_s/(ppa - (1._wp - reps0)*ze_s) 
     745 
    609746   END FUNCTION q_sat_sclr 
    610    !! 
     747 
    611748   FUNCTION q_sat_vctr( pta, ppa,  l_ice ) 
     749 
    612750      REAL(wp), DIMENSION(jpi,jpj) :: q_sat_vctr 
    613751      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta  !: absolute temperature of air [K] 
    … …  
    620758      IF( PRESENT(l_ice) ) lice = l_ice 
    621759      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    622       q_sat_vctr(ji,jj) = q_sat_sclr( pta(ji,jj) , ppa(ji,jj), l_ice=lice ) 
    623       END_2D 
     760         q_sat_vctr(ji,jj) = q_sat_sclr( pta(ji,jj) , ppa(ji,jj), l_ice=lice ) 
     761      END_2D 
     762 
    624763   END FUNCTION q_sat_vctr 
    625    !=============================================================================================== 
    626  
    627  
    628    !=============================================================================================== 
     764 
     765 
    629766   FUNCTION dq_sat_dt_ice_sclr( pta, ppa ) 
    630767      !!--------------------------------------------------------------------------------- 
    … …  
    647784      ! 
    648785   END FUNCTION dq_sat_dt_ice_sclr 
    649    !! 
     786 
    650787   FUNCTION dq_sat_dt_ice_vctr( pta, ppa ) 
     788 
    651789      REAL(wp), DIMENSION(jpi,jpj) :: dq_sat_dt_ice_vctr 
    652790      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta  !: absolute temperature of air [K] 
    … …  
    655793      !!---------------------------------------------------------------------------------- 
    656794      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    657       dq_sat_dt_ice_vctr(ji,jj) = dq_sat_dt_ice_sclr( pta(ji,jj) , ppa(ji,jj) ) 
    658       END_2D 
     795         dq_sat_dt_ice_vctr(ji,jj) = dq_sat_dt_ice_sclr( pta(ji,jj) , ppa(ji,jj) ) 
     796      END_2D 
     797 
    659798   END FUNCTION dq_sat_dt_ice_vctr 
    660    !=============================================================================================== 
    661  
    662  
    663    !=============================================================================================== 
     799 
     800 
    664801   FUNCTION q_air_rh(prha, ptak, ppa) 
    665802      !!---------------------------------------------------------------------------------- 
    … …  
    678815      ! 
    679816      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    680       ze = prha(ji,jj)*e_sat_sclr(ptak(ji,jj)) 
    681       q_air_rh(ji,jj) = ze*reps0/(ppa(ji,jj) - (1. - reps0)*ze) 
     817         ze = prha(ji,jj)*e_sat_sclr(ptak(ji,jj)) 
     818         q_air_rh(ji,jj) = ze*reps0/(ppa(ji,jj) - (1. - reps0)*ze) 
    682819      END_2D 
    683820      ! 
    … …  
    685822 
    686823 
    687    SUBROUTINE UPDATE_QNSOL_TAU( pzu, pTs, pqs, pTa, pqa, pust, ptst, pqst, pwnd, pUb, ppa, prlw, & 
     824   SUBROUTINE UPDATE_QNSOL_TAU( pzu, pTs, pqs, pTa, pqa, pust, ptst, pqst, pwnd, pUb, ppa, prlw, prhoa, & 
    688825      &                         pQns, pTau,  & 
    689826      &                         Qlat) 
    … …  
    705842      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: ppa ! sea-level atmospheric pressure [Pa] 
    706843      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: prlw ! downwelling longwave radiative flux [W/m^2] 
     844      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: prhoa ! air density [kg/m3] 
    707845      ! 
    708846      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQns ! non-solar heat flux to the ocean aka "Qlat + Qsen + Qlw" [W/m^2]] 
    … …  
    715853      !!---------------------------------------------------------------------------------- 
    716854      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    717       zdt = pTa(ji,jj) - pTs(ji,jj) ;  zdt = SIGN( MAX(ABS(zdt),1.E-6_wp), zdt ) 
    718       zdq = pqa(ji,jj) - pqs(ji,jj) ;  zdq = SIGN( MAX(ABS(zdq),1.E-9_wp), zdq ) 
    719       zz0 = pust(ji,jj)/pUb(ji,jj) 
    720       zCd = zz0*zz0 
    721       zCh = zz0*ptst(ji,jj)/zdt 
    722       zCe = zz0*pqst(ji,jj)/zdq 
    723  
    724       CALL BULK_FORMULA( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), zCd, zCh, zCe, & 
    725          &              pwnd(ji,jj), pUb(ji,jj), ppa(ji,jj), & 
    726          &              pTau(ji,jj), zQsen, zQlat ) 
    727  
    728       zQlw = qlw_net_sclr( prlw(ji,jj), pTs(ji,jj) ) ! Net longwave flux 
    729  
    730       pQns(ji,jj) = zQlat + zQsen + zQlw 
    731  
    732       IF( PRESENT(Qlat) ) Qlat(ji,jj) = zQlat 
    733       END_2D 
     855 
     856         zdt = pTa(ji,jj) - pTs(ji,jj) ;  zdt = SIGN( MAX(ABS(zdt),1.E-6_wp), zdt ) 
     857         zdq = pqa(ji,jj) - pqs(ji,jj) ;  zdq = SIGN( MAX(ABS(zdq),1.E-9_wp), zdq ) 
     858         zz0 = pust(ji,jj)/pUb(ji,jj) 
     859         zCd = zz0*zz0 
     860         zCh = zz0*ptst(ji,jj)/zdt 
     861         zCe = zz0*pqst(ji,jj)/zdq 
     862 
     863         CALL BULK_FORMULA( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), zCd, zCh, zCe, & 
     864            &              pwnd(ji,jj), pUb(ji,jj), ppa(ji,jj), prhoa(ji,jj), & 
     865            &              pTau(ji,jj), zQsen, zQlat ) 
     866 
     867         zQlw = qlw_net_sclr( prlw(ji,jj), pTs(ji,jj) ) ! Net longwave flux 
     868 
     869         pQns(ji,jj) = zQlat + zQsen + zQlw 
     870 
     871         IF( PRESENT(Qlat) ) Qlat(ji,jj) = zQlat 
     872 
     873      END_2D 
     874 
    734875   END SUBROUTINE UPDATE_QNSOL_TAU 
    735876 
    … …  
    737878   SUBROUTINE BULK_FORMULA_SCLR( pzu, pTs, pqs, pTa, pqa, & 
    738879      &                          pCd, pCh, pCe,           & 
    739       &                          pwnd, pUb, ppa,         & 
     880      &                          pwnd, pUb, ppa, prhoa,   & 
    740881      &                          pTau, pQsen, pQlat,      & 
    741       &                          pEvap, prhoa, pfact_evap ) 
    742       !!---------------------------------------------------------------------------------- 
    743       REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m) 
    744       REAL(wp), INTENT(in)  :: pTs  ! water temperature at the air-sea interface [K] 
    745       REAL(wp), INTENT(in)  :: pqs  ! satur. spec. hum. at T=pTs   [kg/kg] 
    746       REAL(wp), INTENT(in)  :: pTa  ! potential air temperature at z=pzu [K] 
    747       REAL(wp), INTENT(in)  :: pqa  ! specific humidity at z=pzu [kg/kg] 
     882      &                          pEvap, pfact_evap        ) 
     883      !!---------------------------------------------------------------------------------- 
     884      REAL(wp), INTENT(in)  :: pzu   ! height above the sea-level where all this takes place (normally 10m) 
     885      REAL(wp), INTENT(in)  :: pTs   ! water temperature at the air-sea interface [K] 
     886      REAL(wp), INTENT(in)  :: pqs   ! satur. spec. hum. at T=pTs   [kg/kg] 
     887      REAL(wp), INTENT(in)  :: pTa   ! potential air temperature at z=pzu [K] 
     888      REAL(wp), INTENT(in)  :: pqa   ! specific humidity at z=pzu [kg/kg] 
    748889      REAL(wp), INTENT(in)  :: pCd 
    749890      REAL(wp), INTENT(in)  :: pCh 
    750891      REAL(wp), INTENT(in)  :: pCe 
    751       REAL(wp), INTENT(in)  :: pwnd ! wind speed module at z=pzu [m/s] 
    752       REAL(wp), INTENT(in)  :: pUb  ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s] 
    753       REAL(wp), INTENT(in)  :: ppa ! sea-level atmospheric pressure [Pa] 
     892      REAL(wp), INTENT(in)  :: pwnd  ! wind speed module at z=pzu [m/s] 
     893      REAL(wp), INTENT(in)  :: pUb   ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s] 
     894      REAL(wp), INTENT(in)  :: ppa   ! sea-level atmospheric pressure [Pa] 
     895      REAL(wp), INTENT(in)  :: prhoa ! Air density at z=pzu [kg/m^3] 
    754896      !! 
    755897      REAL(wp), INTENT(out) :: pTau  ! module of the wind stress [N/m^2] 
    … …  
    758900      !! 
    759901      REAL(wp), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s] 
    760       REAL(wp), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3] 
    761902      REAL(wp), INTENT(in) , OPTIONAL :: pfact_evap  ! ABOMINATION: corrective factor for evaporation (doing this against my will! /laurent) 
    762903      !! 
    … …  
    767908      IF( PRESENT(pfact_evap) ) zfact_evap = pfact_evap 
    768909 
    769       !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa") 
    770       ztaa = pTa ! first guess... 
    771       DO jq = 1, 4 
    772          zgamma = gamma_moist( 0.5*(ztaa+pTs) , pqa )  !#LB: why not "0.5*(pqs+pqa)" rather then "pqa" ??? 
    773          ztaa = pTa - zgamma*pzu   ! Absolute temp. is slightly colder... 
    774       END DO 
    775       zrho = rho_air(ztaa, pqa, ppa) 
    776       zrho = rho_air(ztaa, pqa, ppa-zrho*grav*pzu) ! taking into account that we are pzu m above the sea level where SLP is given! 
    777  
    778       zUrho = pUb*MAX(zrho, 1._wp)     ! rho*U10 
     910      zUrho = pUb*MAX(prhoa, 1._wp)     ! rho*U10 
    779911 
    780912      pTau = zUrho * pCd * pwnd ! Wind stress module 
    … …  
    785917 
    786918      IF( PRESENT(pEvap) ) pEvap = - zfact_evap * zevap 
    787       IF( PRESENT(prhoa) ) prhoa = zrho 
    788919 
    789920   END SUBROUTINE BULK_FORMULA_SCLR 
    790    !! 
     921 
    791922   SUBROUTINE BULK_FORMULA_VCTR( pzu, pTs, pqs, pTa, pqa, & 
    792923      &                          pCd, pCh, pCe,           & 
    793       &                          pwnd, pUb, ppa,         & 
     924      &                          pwnd, pUb, ppa, prhoa,   & 
    794925      &                          pTau, pQsen, pQlat,      & 
    795       &                          pEvap, prhoa, pfact_evap ) 
    796       !!---------------------------------------------------------------------------------- 
    797       REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m) 
    798       REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTs  ! water temperature at the air-sea interface [K] 
    799       REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqs  ! satur. spec. hum. at T=pTs   [kg/kg] 
    800       REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTa  ! potential air temperature at z=pzu [K] 
    801       REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqa  ! specific humidity at z=pzu [kg/kg] 
     926      &                          pEvap, pfact_evap ) 
     927      !!---------------------------------------------------------------------------------- 
     928      REAL(wp),                     INTENT(in)  :: pzu   ! height above the sea-level where all this takes place (normally 10m) 
     929      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTs   ! water temperature at the air-sea interface [K] 
     930      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqs   ! satur. spec. hum. at T=pTs   [kg/kg] 
     931      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTa   ! potential air temperature at z=pzu [K] 
     932      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqa   ! specific humidity at z=pzu [kg/kg] 
    802933      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pCd 
    803934      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pCh 
    804935      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pCe 
    805       REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pwnd ! wind speed module at z=pzu [m/s] 
    806       REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pUb  ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s] 
    807       REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: ppa ! sea-level atmospheric pressure [Pa] 
     936      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pwnd  ! wind speed module at z=pzu [m/s] 
     937      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pUb   ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s] 
     938      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: ppa   ! sea-level atmospheric pressure [Pa] 
     939      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: prhoa ! Air density at z=pzu [kg/m^3]  
    808940      !! 
    809941      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pTau  ! module of the wind stress [N/m^2] 
    … …  
    812944      !! 
    813945      REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s] 
    814       REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3] 
    815946      REAL(wp),                     INTENT(in) , OPTIONAL :: pfact_evap  ! ABOMINATION: corrective factor for evaporation (doing this against my will! /laurent) 
    816947      !! 
    … …  
    823954      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    824955 
    825       CALL BULK_FORMULA_SCLR( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), & 
    826          &                    pCd(ji,jj), pCh(ji,jj), pCe(ji,jj),                  & 
    827          &                    pwnd(ji,jj), pUb(ji,jj), ppa(ji,jj),                & 
    828          &                    pTau(ji,jj), pQsen(ji,jj), pQlat(ji,jj),             & 
    829          &                    pEvap=zevap, prhoa=zrho, pfact_evap=zfact_evap       ) 
    830  
    831       IF( PRESENT(pEvap) ) pEvap(ji,jj) = zevap 
    832       IF( PRESENT(prhoa) ) prhoa(ji,jj) = zrho 
    833       END_2D 
     956         CALL BULK_FORMULA_SCLR( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), & 
     957            &                    pCd(ji,jj), pCh(ji,jj), pCe(ji,jj),                  & 
     958            &                    pwnd(ji,jj), pUb(ji,jj), ppa(ji,jj), prhoa(ji,jj),   & 
     959            &                    pTau(ji,jj), pQsen(ji,jj), pQlat(ji,jj),             & 
     960            &                    pEvap=zevap, pfact_evap=zfact_evap                   ) 
     961 
     962         IF( PRESENT(pEvap) ) pEvap(ji,jj) = zevap 
     963      END_2D 
     964 
    834965   END SUBROUTINE BULK_FORMULA_VCTR 
    835966 
    … …  
    847978      !!---------------------------------------------------------------------------------- 
    848979      alpha_sw_vctr = 2.1e-5_wp * MAX(psst(:,:)-rt0 + 3.2_wp, 0._wp)**0.79 
     980 
    849981   END FUNCTION alpha_sw_vctr 
    850982 
    … …  
    861993      !!---------------------------------------------------------------------------------- 
    862994      alpha_sw_sclr = 2.1e-5_wp * MAX(psst-rt0 + 3.2_wp, 0._wp)**0.79 
     995 
    863996   END FUNCTION alpha_sw_sclr 
    864997 
    865998 
    866    !=============================================================================================== 
    867999   FUNCTION qlw_net_sclr( pdwlw, pts,  l_ice ) 
    8681000      !!--------------------------------------------------------------------------------- 
    … …  
    8871019      zt2 = pts*pts 
    8881020      qlw_net_sclr = zemiss*( pdwlw - stefan*zt2*zt2)  ! zemiss used both as the IR albedo and IR emissivity... 
     1021 
    8891022   END FUNCTION qlw_net_sclr 
    890    !! 
     1023 
    8911024   FUNCTION qlw_net_vctr( pdwlw, pts,  l_ice ) 
     1025 
    8921026      REAL(wp), DIMENSION(jpi,jpj) :: qlw_net_vctr 
    8931027      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pdwlw !: downwelling longwave (aka infrared, aka thermal) radiation [W/m^2] 
    … …  
    9001034      IF( PRESENT(l_ice) ) lice = l_ice 
    9011035      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    902       qlw_net_vctr(ji,jj) = qlw_net_sclr( pdwlw(ji,jj) , pts(ji,jj), l_ice=lice ) 
    903       END_2D 
     1036         qlw_net_vctr(ji,jj) = qlw_net_sclr( pdwlw(ji,jj) , pts(ji,jj), l_ice=lice ) 
     1037      END_2D 
     1038 
    9041039   END FUNCTION qlw_net_vctr 
    905    !=============================================================================================== 
     1040 
    9061041 
    9071042   FUNCTION z0_from_Cd( pzu, pCd,  ppsi ) 
     1043 
    9081044      REAL(wp), DIMENSION(jpi,jpj) :: z0_from_Cd        !: roughness length [m] 
    9091045      REAL(wp)                    , INTENT(in) :: pzu   !: reference height zu [m] 
    … …  
    9211057         z0_from_Cd = pzu * EXP( - vkarmn/SQRT(pCd(:,:)) )            !LB: ok, double-checked! 
    9221058      END IF 
     1059 
    9231060   END FUNCTION z0_from_Cd 
    9241061 
     1062 
    9251063   FUNCTION Cd_from_z0( pzu, pz0,  ppsi ) 
     1064 
    9261065      REAL(wp), DIMENSION(jpi,jpj) :: Cd_from_z0        !: (neutral or non-neutral) drag coefficient [] 
    9271066      REAL(wp)                    , INTENT(in) :: pzu   !: reference height zu [m] 
    … …  
    9401079      END IF 
    9411080      Cd_from_z0 = vkarmn2 * Cd_from_z0 * Cd_from_z0 
     1081 
    9421082   END FUNCTION Cd_from_z0 
    9431083 
    … …  
    9651105      ! 
    9661106   END FUNCTION f_m_louis_sclr 
    967    !! 
     1107 
    9681108   FUNCTION f_m_louis_vctr( pzu, pRib, pCdn, pz0 ) 
     1109 
    9691110      REAL(wp), DIMENSION(jpi,jpj)             :: f_m_louis_vctr 
    9701111      REAL(wp),                     INTENT(in) :: pzu 
    … …  
    9731114      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pz0 
    9741115      INTEGER  :: ji, jj 
    975       DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    976       f_m_louis_vctr(ji,jj) = f_m_louis_sclr( pzu, pRib(ji,jj), pCdn(ji,jj), pz0(ji,jj) ) 
    977       END_2D 
     1116 
     1117      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     1118         f_m_louis_vctr(ji,jj) = f_m_louis_sclr( pzu, pRib(ji,jj), pCdn(ji,jj), pz0(ji,jj) ) 
     1119      END_2D 
     1120 
    9781121   END FUNCTION f_m_louis_vctr 
    9791122 
    … …  
    10011144      ! 
    10021145   END FUNCTION f_h_louis_sclr 
    1003    !! 
     1146 
    10041147   FUNCTION f_h_louis_vctr( pzu, pRib, pChn, pz0 ) 
     1148 
    10051149      REAL(wp), DIMENSION(jpi,jpj)             :: f_h_louis_vctr 
    10061150      REAL(wp),                     INTENT(in) :: pzu 
    … …  
    10091153      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pz0 
    10101154      INTEGER  :: ji, jj 
    1011       DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    1012       f_h_louis_vctr(ji,jj) = f_h_louis_sclr( pzu, pRib(ji,jj), pChn(ji,jj), pz0(ji,jj) ) 
    1013       END_2D 
     1155 
     1156      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     1157         f_h_louis_vctr(ji,jj) = f_h_louis_sclr( pzu, pRib(ji,jj), pChn(ji,jj), pz0(ji,jj) ) 
     1158      END_2D 
     1159 
    10141160   END FUNCTION f_h_louis_vctr 
     1161 
    10151162 
    10161163   FUNCTION UN10_from_ustar( pzu, pUzu, pus, ppsi ) 
    … …  
    11231270 
    11241271 
    1125    !!====================================================================== 
     1272 
    11261273END MODULE sbc_phy 

  • NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/OCE/SBC/sbcblk.F90

    r14433 r14592  
    8080   INTEGER , PUBLIC, PARAMETER ::   jp_wndj  =  2   ! index of 10m wind velocity (j-component) (m/s)    at T-point 
    8181   INTEGER , PUBLIC, PARAMETER ::   jp_tair  =  3   ! index of 10m air temperature             (Kelvin) 
    82    INTEGER , PUBLIC, PARAMETER ::   jp_humi  =  4   ! index of specific humidity               ( % ) 
     82   INTEGER , PUBLIC, PARAMETER ::   jp_humi  =  4   ! index of specific humidity               (kg/kg) 
    8383   INTEGER , PUBLIC, PARAMETER ::   jp_qsr   =  5   ! index of solar heat                      (W/m2) 
    8484   INTEGER , PUBLIC, PARAMETER ::   jp_qlw   =  6   ! index of Long wave                       (W/m2) 
    … …  
    503503      INTEGER, INTENT(in) ::   kt   ! ocean time step 
    504504      !!---------------------------------------------------------------------- 
    505       REAL(wp), DIMENSION(jpi,jpj) ::   zssq, zcd_du, zsen, zlat, zevp 
     505      REAL(wp), DIMENSION(jpi,jpj) ::   zssq, zcd_du, zsen, zlat, zevp, zpre, ztheta 
    506506      REAL(wp) :: ztst 
    507507      LOGICAL  :: llerr 
    … …  
    540540      IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 
    541541 
    542          ! Specific humidity of air at z=rn_zqt ! 
     542         ! Specific humidity of air at z=rn_zqt 
    543543         SELECT CASE( nhumi ) 
    544544         CASE( np_humi_sph ) 
    … …  
    553553         END SELECT 
    554554 
    555          ! POTENTIAL temperature of air at z=rn_zqt 
    556          !      based on adiabatic lapse-rate (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2 
    557          !      (most reanalysis products provide absolute temp., not potential temp.) 
     555         ! Potential temperature of air at z=rn_zqt (most reanalysis products provide absolute temp., not potential temp.) 
    558556         IF( ln_tair_pot ) THEN 
    559             ! temperature read into file is already potential temperature, do nothing... 
    560             theta_air_zt(:,:) = sf(jp_tair )%fnow(:,:,1) 
     557            theta_air_zt(:,:)        = sf(jp_tair )%fnow(:,:,1) 
    561558         ELSE 
    562559            ! temperature read into file is ABSOLUTE temperature (that's the case for ECMWF products for example...) 
    563560            IF((kt==nit000).AND.lwp) WRITE(numout,*) ' *** sbc_blk() => air temperature converted from ABSOLUTE to POTENTIAL!' 
    564             theta_air_zt(:,:) = sf(jp_tair )%fnow(:,:,1) + gamma_moist( sf(jp_tair )%fnow(:,:,1), q_air_zt(:,:) ) * rn_zqt 
     561            zpre(:,:)         = pres_temp( q_air_zt(:,:), sf(jp_slp)%fnow(:,:,1), rn_zu, pta=sf(jp_tair)%fnow(:,:,1) ) 
     562            theta_air_zt(:,:) = theta_exner( sf(jp_tair)%fnow(:,:,1), zpre(:,:) ) 
    565563         ENDIF 
    566564         ! 
    … …  
    588586         tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1)    !#LB: should it be POTENTIAL temperature instead ???? 
    589587         !tatm_ice(:,:) = theta_air_zt(:,:)         !#LB: THIS! ? 
    590  
    591          qatm_ice(:,:) = q_air_zt(:,:) !#LB: 
    592  
    593          tprecip(:,:)     = sf(jp_prec)%fnow(:,:,1) * rn_pfac 
    594          sprecip(:,:)     = sf(jp_snow)%fnow(:,:,1) * rn_pfac 
    595          wndi_ice(:,:)    = sf(jp_wndi)%fnow(:,:,1) 
    596          wndj_ice(:,:)    = sf(jp_wndj)%fnow(:,:,1) 
     588         qatm_ice(:,:) = q_air_zt(:,:) 
     589         tprecip(:,:)  = sf(jp_prec)%fnow(:,:,1) * rn_pfac 
     590         sprecip(:,:)  = sf(jp_snow)%fnow(:,:,1) * rn_pfac 
     591         wndi_ice(:,:) = sf(jp_wndi)%fnow(:,:,1) 
     592         wndj_ice(:,:) = sf(jp_wndj)%fnow(:,:,1) 
    597593      ENDIF 
    598594#endif 
    … …  
    652648      REAL(wp), DIMENSION(jpi,jpj) ::   zch_oce           ! sensible heat transfert coefficient over ocean 
    653649      REAL(wp), DIMENSION(jpi,jpj) ::   zce_oce           ! latent   heat transfert coefficient over ocean 
     650      REAL(wp), DIMENSION(jpi,jpj) ::   zsspt             ! potential sea-surface temperature [K] 
     651      REAL(wp), DIMENSION(jpi,jpj) ::   zpre, ztabs 
    654652      REAL(wp), DIMENSION(jpi,jpj) ::   zztmp1, zztmp2 
    655653      !!--------------------------------------------------------------------- 
    … …  
    658656      !                           ! Temporary conversion from Celcius to Kelvin (and set minimum value far above 0 K) 
    659657      ptsk(:,:) = pst(:,:) + rt0  ! by default: skin temperature = "bulk SST" (will remain this way if NCAR algorithm used!) 
     658 
     659      ! sea surface potential temperature [K] 
     660      zsspt(:,:) = theta_exner( ptsk(:,:), pslp(:,:) ) 
    660661 
    661662      ! --- cloud cover --- ! 
    … …  
    705706      IF( ln_skin_cs .OR. ln_skin_wl ) THEN 
    706707         !! Backup "bulk SST" and associated spec. hum. 
    707          zztmp1(:,:) = ptsk(:,:) 
     708         zztmp1(:,:) = zsspt(:,:) 
    708709         zztmp2(:,:) = pssq(:,:) 
    709710      ENDIF 
    … …  
    714715 
    715716      CASE( np_NCAR      ) 
    716          CALL turb_ncar    (     rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, & 
     717         CALL turb_ncar    (     rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, & 
    717718            &                zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu , & 
    718719            &                nb_iter=nn_iter_algo ) 
    719720         ! 
    720721      CASE( np_COARE_3p0 ) 
    721          CALL turb_coare3p0( kt, rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, & 
     722         CALL turb_coare3p0( kt, rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, & 
    722723            &                ln_skin_cs, ln_skin_wl,                            & 
    723724            &                zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu,  & 
    … …  
    726727         ! 
    727728      CASE( np_COARE_3p6 ) 
    728          CALL turb_coare3p6( kt, rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, & 
     729         CALL turb_coare3p6( kt, rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, & 
    729730            &                ln_skin_cs, ln_skin_wl,                            & 
    730731            &                zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu,  & 
    … …  
    733734         ! 
    734735      CASE( np_ECMWF     ) 
    735          CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, & 
     736         CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, & 
    736737            &                ln_skin_cs, ln_skin_wl,                            & 
    737738            &                zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu,  & 
    … …  
    740741         ! 
    741742      CASE( np_ANDREAS   ) 
    742          CALL turb_andreas (     rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, & 
     743         CALL turb_andreas (     rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, & 
    743744            &                zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu , & 
    744745            &                nb_iter=nn_iter_algo   ) 
    … …  
    762763      IF( ln_skin_cs .OR. ln_skin_wl ) THEN 
    763764         !! ptsk and pssq have been updated!!! 
     765         ptsk(:,:) = ptsk(:,:) + ( zsspt(:,:) - zztmp1(:,:) ) 
    764766         !! 
    765767         !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of ptsk and pssq: 
    766768         WHERE ( fr_i(:,:) > 0.001_wp ) 
    767769            ! sea-ice present, we forget about the update, using what we backed up before call to turb_*() 
    768             ptsk(:,:) = zztmp1(:,:) 
    769             pssq(:,:) = zztmp2(:,:) 
     770            zsspt(:,:) = zztmp1(:,:) 
     771            pssq(:,:)  = zztmp2(:,:) 
    770772         END WHERE 
    771773      END IF 
    … …  
    775777 
    776778      IF( ln_abl ) THEN         !==  ABL formulation  ==!   multiplication by rho_air and turbulent fluxes computation done in ablstp 
     779 
    777780         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    778781            zztmp = zU_zu(ji,jj) 
    … …  
    781784            psen(ji,jj)   = zztmp * zch_oce(ji,jj) 
    782785            pevp(ji,jj)   = zztmp * zce_oce(ji,jj) 
    783             rhoa(ji,jj)   = rho_air( ptair(ji,jj), pqair(ji,jj), pslp(ji,jj) ) 
     786            zpre(ji,jj)   = pres_temp( pqair(ji,jj), pslp(ji,jj), rn_zu, ptpot=ptair(ji,jj), pta=ztabs(ji,jj) ) 
     787            rhoa(ji,jj)   = rho_air( ztabs(ji,jj), pqair(ji,jj), zpre(ji,jj) ) 
    784788         END_2D 
     789 
    785790      ELSE                      !==  BLK formulation  ==!   turbulent fluxes computation 
    786          CALL BULK_FORMULA( rn_zu, ptsk(:,:), pssq(:,:), theta_zu(:,:), q_zu(:,:), & 
     791 
     792         zpre(:,:) = pres_temp( q_zu(:,:), pslp(:,:), rn_zu, ptpot=theta_zu(:,:), pta=ztabs(:,:) ) 
     793         rhoa(:,:) = rho_air( ztabs(:,:), q_zu(:,:), zpre(:,:) ) 
     794         !!GS: for debug, to be removed  
     795         CALL iom_put( "pres_zu", zpre(:,:)*tmask(:,:,1) ) 
     796         CALL iom_put( "slp", pslp(:,:)*tmask(:,:,1) ) 
     797         CALL iom_put( "tabs_zu", (ztabs(:,:)-rt0)*tmask(:,:,1) ) 
     798 
     799         CALL BULK_FORMULA( rn_zu, zsspt(:,:), pssq(:,:), theta_zu(:,:), q_zu(:,:), & 
    787800            &               zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:),          & 
    788             &               wndm(:,:), zU_zu(:,:), pslp(:,:),                  & 
     801            &               wndm(:,:), zU_zu(:,:), pslp(:,:), rhoa(:,:),       & 
    789802            &               taum(:,:), psen(:,:), plat(:,:),                   & 
    790             &               pEvap=pevp(:,:), prhoa=rhoa(:,:), pfact_evap=rn_efac ) 
     803            &               pEvap=pevp(:,:), pfact_evap=rn_efac ) 
    791804 
    792805         psen(:,:) = psen(:,:) * tmask(:,:,1) 
    … …  
    851864         ENDIF 
    852865         ! 
    853       ENDIF !IF( ln_abl ) 
     866      ENDIF ! ln_blk / ln_abl 
    854867 
    855868      ptsk(:,:) = ( ptsk(:,:) - rt0 ) * tmask(:,:,1)  ! Back to Celsius 
    856869 
     870      CALL iom_put( "sspt", zsspt(:,:)-rt0 ) 
    857871      IF( ln_skin_cs .OR. ln_skin_wl ) THEN 
    858872         CALL iom_put( "t_skin" ,  ptsk        )  ! T_skin in Celsius 
    … …  
    970984      !! ** Method  :   compute momentum using bulk formulation 
    971985      !!                formulea, ice variables and read atmospheric fields. 
    972       !!                NB: ice drag coefficient is assumed to be a constant 
    973986      !!--------------------------------------------------------------------- 
    974987      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pslp    ! sea-level pressure [Pa] 
    975988      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pwndi   ! atmospheric wind at T-point [m/s] 
    976989      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pwndj   ! atmospheric wind at T-point [m/s] 
    977       REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   ptair   ! atmospheric wind at T-point [m/s] 
    978       REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pqair   ! atmospheric wind at T-point [m/s] 
     990      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   ptair   ! atmospheric potential temperature at T-point [K] 
     991      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pqair   ! atmospheric specific humidity at T-point [kg/kg] 
    979992      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   puice   ! sea-ice velocity on I or C grid [m/s] 
    980993      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pvice   ! " 
    … …  
    9901003      REAL(wp) ::   zootm_su                      ! sea-ice surface mean temperature 
    9911004      REAL(wp) ::   zztmp1, zztmp2                ! temporary scalars 
    992       REAL(wp), DIMENSION(jpi,jpj) :: ztmp        ! temporary array 
    993       !!--------------------------------------------------------------------- 
    994       ! 
    995       ! LB: ptsui is in K !!! 
     1005      REAL(wp), DIMENSION(jpi,jpj) :: ztmp, zsipt ! temporary array 
     1006      !!--------------------------------------------------------------------- 
    9961007      ! 
    9971008      ! ------------------------------------------------------------ ! 
    … …  
    10001011      ! C-grid ice dynamics :   U & V-points (same as ocean) 
    10011012      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    1002       wndm_ice(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) ) 
     1013         wndm_ice(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) ) 
    10031014      END_2D 
    10041015      ! 
    1005       ! Make ice-atm. drag dependent on ice concentration 
    1006  
    1007  
     1016      ! potential sea-ice surface temperature [K] 
     1017      zsipt(:,:) = theta_exner( ptsui(:,:), pslp(:,:) ) 
     1018       
     1019      ! sea-ice <-> atmosphere bulk transfer coefficients 
    10081020      SELECT CASE( nblk_ice ) 
    10091021 
    … …  
    10191031      CASE( np_ice_an05 )  ! calculate new drag from Lupkes(2015) equations 
    10201032         ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ 
    1021          CALL turb_ice_an05( rn_zqt, rn_zu, ptsui, ptair, ztmp, pqair, wndm_ice,       & 
     1033         CALL turb_ice_an05( rn_zqt, rn_zu, zsipt, ptair, ztmp, pqair, wndm_ice,       & 
    10221034            &                      Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i ) 
    10231035         !! 
    10241036      CASE( np_ice_lu12 ) 
    10251037         ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ 
    1026          CALL turb_ice_lu12( rn_zqt, rn_zu, ptsui, ptair, ztmp, pqair, wndm_ice, fr_i, & 
     1038         CALL turb_ice_lu12( rn_zqt, rn_zu, zsipt, ptair, ztmp, pqair, wndm_ice, fr_i, & 
    10271039            &                      Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i ) 
    10281040         !! 
    10291041      CASE( np_ice_lg15 )  ! calculate new drag from Lupkes(2015) equations 
    10301042         ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ 
    1031          CALL turb_ice_lg15( rn_zqt, rn_zu, ptsui, ptair, ztmp, pqair, wndm_ice, fr_i, & 
     1043         CALL turb_ice_lg15( rn_zqt, rn_zu, zsipt, ptair, ztmp, pqair, wndm_ice, fr_i, & 
    10321044            &                      Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i ) 
    10331045         !! 
    10341046      END SELECT 
    10351047 
    1036       IF( iom_use('Cd_ice').OR.iom_use('Ce_ice').OR.iom_use('Ch_ice').OR.iom_use('taum_ice').OR.iom_use('utau_ice').OR.iom_use('vtau_ice') ) & 
    1037          & ztmp(:,:) = ( 1._wp - MAX(0._wp, SIGN( 1._wp, 1.E-6_wp - fr_i )) )*tmask(:,:,1) ! mask for presence of ice ! 
    1038  
     1048      ztmp(:,:) = ( 1._wp - MAX(0._wp, SIGN( 1._wp, 1.E-6_wp - fr_i )) )*tmask(:,:,1) ! mask for presence of ice ! 
     1049      IF( iom_use('spt_ice')) CALL iom_put("spt_ice", (zsipt-rt0)*ztmp) 
    10391050      IF( iom_use('Cd_ice') ) CALL iom_put("Cd_ice", Cd_ice*ztmp) 
    10401051      IF( iom_use('Ce_ice') ) CALL iom_put("Ce_ice", Ce_ice*ztmp) 
    … …  
    10491060         DO_2D( 0, 1, 0, 1 )    ! at T point 
    10501061            zztmp1        = rhoa(ji,jj) * Cd_ice(ji,jj) * wndm_ice(ji,jj) 
    1051             putaui(ji,jj) =  zztmp1 * pwndi(ji,jj) 
    1052             pvtaui(ji,jj) =  zztmp1 * pwndj(ji,jj) 
     1062            putaui(ji,jj) = zztmp1 * pwndi(ji,jj) 
     1063            pvtaui(ji,jj) = zztmp1 * pwndj(ji,jj) 
    10531064         END_2D 
    10541065 
    … …  
    10731084            &                               , tab2d_2=pvtaui  , clinfo2='          pvtaui : ' ) 
    10741085      ELSE ! ln_abl 
     1086 
    10751087         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    1076          pcd_dui(ji,jj) = wndm_ice(ji,jj) * Cd_ice(ji,jj) 
    1077          pseni  (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj) 
    1078          pevpi  (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj) 
     1088            pcd_dui(ji,jj) = wndm_ice(ji,jj) * Cd_ice(ji,jj) 
     1089            pseni  (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj) 
     1090            pevpi  (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj) 
    10791091         END_2D 
    1080          !#LB: 
    10811092         pssqi(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ; ! more accurate way to obtain ssq ! 
    1082          !#LB. 
    1083       ENDIF !IF( ln_blk ) 
     1093 
     1094      ENDIF ! ln_blk  / ln_abl 
    10841095      ! 
    10851096      IF(sn_cfctl%l_prtctl)  CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice: wndm_ice : ') 
    … …  
    11121123      !! 
    11131124      INTEGER  ::   ji, jj, jl               ! dummy loop indices 
    1114       REAL(wp) ::   zst, zst3, zsq           ! local variable 
     1125      REAL(wp) ::   zst, zst3, zsq, zsipt    ! local variable 
    11151126      REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      - 
    11161127      REAL(wp) ::   zztmp, zzblk, zztmp1, z1_rLsub   !   -      - 
    … …  
    11261137      zcoef_dqlw = 4._wp * emiss_i * stefan             ! local scalars 
    11271138      ! 
    1128  
    11291139      zztmp = 1. / ( 1. - albo ) 
    11301140      dqla_ice(:,:,:) = 0._wp 
    1131  
    11321141      !                                     ! ========================== ! 
    11331142      DO jl = 1, jpl                        !  Loop over ice categories  ! 
    … …  
    11351144         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    11361145 
    1137                zst = ptsu(ji,jj,jl)                           ! surface temperature of sea-ice [K] 
    1138                zsq = q_sat( zst, pslp(ji,jj), l_ice=.TRUE. )  ! surface saturation specific humidity when ice present 
    1139  
    1140                ! ----------------------------! 
    1141                !      I   Radiative FLUXES   ! 
    1142                ! ----------------------------! 
    1143                ! Short Wave (sw) 
    1144                qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj) 
    1145  
    1146                ! Long  Wave (lw) 
    1147                zst3 = zst * zst * zst 
    1148                z_qlw(ji,jj,jl)   = emiss_i * ( pdqlw(ji,jj) - stefan * zst * zst3 ) * tmask(ji,jj,1) 
    1149                ! lw sensitivity 
    1150                z_dqlw(ji,jj,jl)  = zcoef_dqlw * zst3 
    1151  
    1152                ! ----------------------------! 
    1153                !     II    Turbulent FLUXES  ! 
    1154                ! ----------------------------! 
    1155  
    1156                ! ... turbulent heat fluxes with Ch_ice recalculated in blk_ice_1 
    1157  
    1158                ! Common term in bulk F. equations... 
    1159                zzblk = rhoa(ji,jj) * wndm_ice(ji,jj) 
    1160  
    1161                ! Sensible Heat 
    1162                zztmp1 = zzblk * rCp_air * Ch_ice(ji,jj) 
    1163                z_qsb (ji,jj,jl) = zztmp1 * (zst - theta_zu_i(ji,jj)) 
    1164                z_dqsb(ji,jj,jl) = zztmp1                        ! ==> Qsens sensitivity (Dqsb_ice/Dtn_ice) 
    1165  
    1166                ! Latent Heat 
    1167                zztmp1 = zzblk * rLsub * Ce_ice(ji,jj) 
    1168                qla_ice(ji,jj,jl) = MAX( zztmp1 * (zsq - q_zu_i(ji,jj)) , 0._wp )   ! #LB: only sublimation (and not condensation) ??? 
    1169                IF(qla_ice(ji,jj,jl)>0._wp) dqla_ice(ji,jj,jl) = zztmp1*dq_sat_dt_ice(zst, pslp(ji,jj)) ! ==> Qlat sensitivity  (dQlat/dT) 
    1170                !                                                                                       !#LB: dq_sat_dt_ice() in "sbc_phy.F90" 
    1171                !#LB: without this unjustified "condensation sensure": 
    1172                !qla_ice( ji,jj,jl) = zztmp1 * (zsq - q_zu_i(ji,jj)) 
    1173                !dqla_ice(ji,jj,jl) = zztmp1 * dq_sat_dt_ice(zst, pslp(ji,jj)) ! ==> Qlat sensitivity  (dQlat/dT) 
    1174  
    1175  
    1176                ! ----------------------------! 
    1177                !     III    Total FLUXES     ! 
    1178                ! ----------------------------! 
    1179                ! Downward Non Solar flux 
    1180                qns_ice (ji,jj,jl) =     z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl) 
    1181                ! Total non solar heat flux sensitivity for ice 
    1182                dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) ) !#LB: correct signs ???? 
     1146            zst   = ptsu(ji,jj,jl)                                ! surface temperature of sea-ice [K] 
     1147            zsq   = q_sat( zst, pslp(ji,jj), l_ice=.TRUE. )       ! surface saturation specific humidity when ice present 
     1148            zsipt = theta_exner( zst, pslp(ji,jj) )               ! potential sea-ice surface temperature [K]   
     1149 
     1150            ! ----------------------------! 
     1151            !      I   Radiative FLUXES   ! 
     1152            ! ----------------------------! 
     1153            ! Short Wave (sw) 
     1154            qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj) 
     1155 
     1156            ! Long  Wave (lw) 
     1157            zst3 = zst * zst * zst 
     1158            z_qlw(ji,jj,jl)   = emiss_i * ( pdqlw(ji,jj) - stefan * zst * zst3 ) * tmask(ji,jj,1) 
     1159            ! lw sensitivity 
     1160            z_dqlw(ji,jj,jl)  = zcoef_dqlw * zst3 
     1161 
     1162            ! ----------------------------! 
     1163            !     II    Turbulent FLUXES  ! 
     1164            ! ----------------------------! 
     1165 
     1166            ! ... turbulent heat fluxes with Ch_ice recalculated in blk_ice_1 
     1167 
     1168            ! Common term in bulk F. equations... 
     1169            zzblk = rhoa(ji,jj) * wndm_ice(ji,jj) 
     1170 
     1171            ! Sensible Heat 
     1172            zztmp1 = zzblk * rCp_air * Ch_ice(ji,jj) 
     1173            z_qsb (ji,jj,jl) = zztmp1 * (zsipt - theta_zu_i(ji,jj)) 
     1174            z_dqsb(ji,jj,jl) = zztmp1                        ! ==> Qsens sensitivity (Dqsb_ice/Dtn_ice) 
     1175 
     1176            ! Latent Heat 
     1177            zztmp1 = zzblk * rLsub * Ce_ice(ji,jj) 
     1178            qla_ice(ji,jj,jl) = MAX( zztmp1 * (zsq - q_zu_i(ji,jj)) , 0._wp )   ! #LB: only sublimation (and not condensation) ??? 
     1179            IF(qla_ice(ji,jj,jl)>0._wp) dqla_ice(ji,jj,jl) = zztmp1*dq_sat_dt_ice(zst, pslp(ji,jj)) ! ==> Qlat sensitivity  (dQlat/dT) 
     1180            !                                                                                       !#LB: dq_sat_dt_ice() in "sbc_phy.F90" 
     1181            !#LB: without this unjustified "condensation sensure": 
     1182            !qla_ice( ji,jj,jl) = zztmp1 * (zsq - q_zu_i(ji,jj)) 
     1183            !dqla_ice(ji,jj,jl) = zztmp1 * dq_sat_dt_ice(zst, pslp(ji,jj)) ! ==> Qlat sensitivity  (dQlat/dT) 
     1184 
     1185 
     1186            ! ----------------------------! 
     1187            !     III    Total FLUXES     ! 
     1188            ! ----------------------------! 
     1189            ! Downward Non Solar flux 
     1190            qns_ice (ji,jj,jl) =     z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl) 
     1191            ! Total non solar heat flux sensitivity for ice 
     1192            dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) ) !#LB: correct signs ???? 
    11831193 
    11841194         END_2D 

  • NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/OCE/SBC/sbcblk_algo_coare3p0.F90

    r14072 r14592  
    189189      REAL(wp), DIMENSION(jpi,jpj) :: zeta_u        ! stability parameter at height zu 
    190190      REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2 
     191      REAL(wp), DIMENSION(jpi,jpj) :: zpre, zrhoa, zta 
    191192      ! 
    192193      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zeta_t  ! stability parameter at height zt 
    … …  
    321322         ENDIF 
    322323 
     324         IF(( l_use_cs ).OR.( l_use_wl )) THEN 
     325            zpre(:,:) = pres_temp( q_zu(:,:), slp(:,:), zu, ptpot=t_zu(:,:), pta=zta(:,:) ) 
     326            zrhoa(:,:) = rho_air( zta(:,:), q_zu(:,:), zpre(:,:) ) 
     327         ENDIF 
     328 
    323329         IF( l_use_cs ) THEN 
    324330            !! Cool-skin contribution 
    325331 
    326             CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, & 
    327                &                   ztmp1, zeta_u,  Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> zeta_u 
     332            CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, zrhoa, & 
     333               &                   ztmp1, zeta_u, Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> zeta_u 
    328334 
    329335            CALL CS_COARE( Qsw, ztmp1, u_star, zsst, ztmp2 )  ! ! Qnsol -> ztmp1 / Qlat -> ztmp2 
    … …  
    336342         IF( l_use_wl ) THEN 
    337343            !! Warm-layer contribution 
    338             CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, & 
     344            CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, zrhoa, & 
    339345               &                   ztmp1, zeta_u)  ! Qnsol -> ztmp1 / Tau -> zeta_u 
    340346            !! In WL_COARE or , Tau_ac and Qnt_ac must be updated at the final itteration step => add a flag to do this! 

  • NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/OCE/SBC/sbcblk_algo_coare3p6.F90

    r14072 r14592  
    179179      REAL(wp), DIMENSION(jpi,jpj) :: zeta_u        ! stability parameter at height zu 
    180180      REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2 
     181      REAL(wp), DIMENSION(jpi,jpj) :: zpre, zrhoa, zta 
    181182      ! 
    182183      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zeta_t  ! stability parameter at height zt 
    … …  
    311312         ENDIF 
    312313 
     314         IF(( l_use_cs ).OR.( l_use_wl )) THEN 
     315            zpre(:,:)  = pres_temp( q_zu(:,:), slp(:,:), zu, ptpot=t_zu(:,:), pta=zta(:,:) ) 
     316            zrhoa(:,:) = rho_air( zta(:,:), q_zu(:,:), zpre(:,:) ) 
     317         ENDIF 
     318 
    313319         IF( l_use_cs ) THEN 
    314320            !! Cool-skin contribution 
    315321 
    316             CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, & 
    317                &                   ztmp1, zeta_u,  Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> zeta_u 
     322            CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, zrhoa, & 
     323               &                   ztmp1, zeta_u, Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> zeta_u 
    318324 
    319325            CALL CS_COARE( Qsw, ztmp1, u_star, zsst, ztmp2 )  ! ! Qnsol -> ztmp1 / Qlat -> ztmp2 
    … …  
    326332         IF( l_use_wl ) THEN 
    327333            !! Warm-layer contribution 
    328             CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, & 
     334            CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, zrhoa, & 
    329335               &                   ztmp1, zeta_u)  ! Qnsol -> ztmp1 / Tau -> zeta_u 
    330336            !! In WL_COARE or , Tau_ac and Qnt_ac must be updated at the final itteration step => add a flag to do this! 

  • NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/OCE/SBC/sbcblk_algo_ecmwf.F90

    r14072 r14592  
    182182      REAL(wp), DIMENSION(jpi,jpj) :: Linv  !: 1/L (inverse of Monin Obukhov length... 
    183183      REAL(wp), DIMENSION(jpi,jpj) :: z0, z0t, z0q 
     184      REAL(wp), DIMENSION(jpi,jpj) :: zrhoa, zpre, zta 
    184185      ! 
    185186      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zsst  ! to back up the initial bulk SST 
    … …  
    336337         func_h = log(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv) 
    337338 
     339         IF(( l_use_cs ).OR.( l_use_wl )) THEN 
     340            zpre(:,:)  = pres_temp( q_zu(:,:), slp(:,:), zu, ptpot=t_zu(:,:), pta=zta(:,:) ) 
     341            zrhoa(:,:) = rho_air( zta(:,:), q_zu(:,:), zpre(:,:) ) 
     342         ENDIF 
    338343 
    339344         IF( l_use_cs ) THEN 
    340345            !! Cool-skin contribution 
    341346 
    342             CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, & 
    343                &                   ztmp1, ztmp0,  Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> ztmp0 
     347 
     348            CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, zrhoa, & 
     349               &                   ztmp1, ztmp0, Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> ztmp0 
    344350 
    345351            CALL CS_ECMWF( Qsw, ztmp1, u_star, zsst )  ! Qnsol -> ztmp1 
    … …  
    353359         IF( l_use_wl ) THEN 
    354360            !! Warm-layer contribution 
    355             CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, & 
     361            CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, zrhoa, & 
    356362               &                   ztmp1, ztmp2)  ! Qnsol -> ztmp1 / Tau -> ztmp2 
    357363            CALL WL_ECMWF( Qsw, ztmp1, u_star, zsst ) 
Note: See TracChangeset for help on using the changeset viewer.