119        undef => const_undef, &
 
  120        pre00 => const_pre00, &
 
  121        tem00 => const_tem00, &
 
  122        rdry  => const_rdry,  &
 
  123        cpdry => const_cpdry, &
 
  124        rvap  => const_rvap,  &
 
  127        qsat => atmos_saturation_dens2qsat_all
 
  130        atmos_hydrometeor_lhv, &
 
  131        atmos_hydrometeor_lhs, &
 
  137        bulkflux_diagnose_surface
 
  140     integer,          
intent(in)    :: IA, IS, IE
 
  141     integer,          
intent(in)    :: JA, JS, JE
 
  142     real(RP),         
intent(in)    :: TMPA     (IA,JA)                     
 
  143     real(RP),         
intent(in)    :: PRSA     (IA,JA)                     
 
  144     real(RP),         
intent(in)    :: WA       (IA,JA)                     
 
  145     real(RP),         
intent(in)    :: UA       (IA,JA)                     
 
  146     real(RP),         
intent(in)    :: VA       (IA,JA)                     
 
  147     real(RP),         
intent(in)    :: RHOA     (IA,JA)                     
 
  148     real(RP),         
intent(in)    :: QVA      (IA,JA)                     
 
  149     real(RP),         
intent(in)    :: LH       (IA,JA)                     
 
  150     real(RP),         
intent(in)    :: Z1       (IA,JA)                     
 
  151     real(RP),         
intent(in)    :: PBL      (IA,JA)                     
 
  152     real(RP),         
intent(in)    :: RHOS     (IA,JA)                     
 
  153     real(RP),         
intent(in)    :: PRSS     (IA,JA)                     
 
  154     real(RP),         
intent(in)    :: RFLXD    (IA,JA,N_RAD_DIR,N_RAD_RGN) 
 
  155     real(RP),         
intent(in)    :: TG       (IA,JA)                     
 
  156     real(RP),         
intent(in)    :: WSTR     (IA,JA)                     
 
  157     real(RP),         
intent(in)    :: QVEF     (IA,JA)                     
 
  158     real(RP),         
intent(in)    :: ALBEDO   (IA,JA,N_RAD_DIR,N_RAD_RGN) 
 
  159     real(RP),         
intent(in)    :: Rb       (IA,JA)                     
 
  160     real(RP),         
intent(in)    :: TC_dZ    (IA,JA)                     
 
  161     real(RP),         
intent(in)    :: Z0M      (IA,JA)                     
 
  162     real(RP),         
intent(in)    :: Z0H      (IA,JA)                     
 
  163     real(RP),         
intent(in)    :: Z0E      (IA,JA)                     
 
  164     logical,          
intent(in)    :: calc_flag(IA,JA)                     
 
  165     real(DP),         
intent(in)    :: dt                                   
 
  166     character(len=*), 
intent(in)    :: model_name
 
  168     real(RP),         
intent(inout) :: TMPS     (IA,JA)                     
 
  170     real(RP),         
intent(out)   :: ZMFLX    (IA,JA)                     
 
  171     real(RP),         
intent(out)   :: XMFLX    (IA,JA)                     
 
  172     real(RP),         
intent(out)   :: YMFLX    (IA,JA)                     
 
  173     real(RP),         
intent(out)   :: SHFLX    (IA,JA)                     
 
  174     real(RP),         
intent(out)   :: LHFLX    (IA,JA)                     
 
  175     real(RP),         
intent(out)   :: QVFLX    (IA,JA)                     
 
  176     real(RP),         
intent(out)   :: GFLX     (IA,JA)                     
 
  177     real(RP),         
intent(out)   :: Ustar    (IA,JA)                     
 
  178     real(RP),         
intent(out)   :: Tstar    (IA,JA)                     
 
  179     real(RP),         
intent(out)   :: Qstar    (IA,JA)                     
 
  180     real(RP),         
intent(out)   :: Wstar    (IA,JA)                     
 
  181     real(RP),         
intent(out)   :: RLmo     (IA,JA)                     
 
  182     real(RP),         
intent(out)   :: U10      (IA,JA)                     
 
  183     real(RP),         
intent(out)   :: V10      (IA,JA)                     
 
  184     real(RP),         
intent(out)   :: T2       (IA,JA)                     
 
  185     real(RP),         
intent(out)   :: Q2       (IA,JA)                     
 
  187     real(RP), 
parameter :: dTS0     = 1.0e-4_rp 
 
  188     real(RP), 
parameter :: redf_min = 1.0e-2_rp 
 
  189     real(RP), 
parameter :: redf_max = 1.0e+0_rp 
 
  190     real(RP), 
parameter :: TFa      = 0.5e+0_rp 
 
  191     real(RP), 
parameter :: TFb      = 1.1e+0_rp 
 
  193     real(RP) :: TMPS1(IA,JA)
 
  214     real(RP) :: Uabs, dUabs 
 
  217     real(RP) :: QVsat, dQVsat    
 
  218     real(RP) :: QVS(IA,JA), dQVS 
 
  222     real(RP) :: FracU10(IA,JA), dFracU10 
 
  223     real(RP) :: FracT2 (IA,JA), dFracT2  
 
  224     real(RP) :: FracQ2 (IA,JA), dFracQ2  
 
  231     log_progress(*) 
'coupler / physics / surface / SKIN' 
  237        tmps1(i,j) = tmps(i,j)
 
  260        if ( calc_flag(i,j) ) 
then 
  266           oldres = huge(0.0_rp)
 
  267           olddts = cpl_phy_sfc_skin_dts_max * dt
 
  270           do n = 1, cpl_phy_sfc_skin_itr_max
 
  272              call qsat( tmps1(i,j),      rhos(i,j), qvsat  )
 
  273              call qsat( tmps1(i,j)+dts0, rhos(i,j), dqvsat )
 
  277              qvs(i,j) = ( 1.0_rp-qvef(i,j) ) * qva(i,j) &
 
  278                       + (        qvef(i,j) ) * qvsat
 
  279              dqvs     = ( 1.0_rp-qvef(i,j) ) * qva(i,j) &
 
  280                       + (        qvef(i,j) ) * dqvsat
 
  282              uabs = sqrt( wa(i,j)**2 + ua(i,j)**2 + va(i,j)**2 )
 
  284              call bulkflux( tmpa(i,j), tmps1(i,j),                 & 
 
  285                             prsa(i,j), prss(i,j),                 & 
 
  286                             qva(i,j), qvs(i,j),                 & 
 
  287                             uabs, z1(i,j), pbl(i,j),               & 
 
  288                             z0m(i,j), z0h(i,j), z0e(i,j),          & 
 
  289                             ustar(i,j), tstar(i,j), qstar(i,j),    & 
 
  290                             wstar(i,j), rlmo(i,j), ra,             & 
 
  291                             fracu10(i,j), fract2(i,j), fracq2(i,j) ) 
 
  293              call bulkflux( tmpa(i,j), tmps1(i,j)+dts0,  & 
 
  294                             prsa(i,j), prss(i,j),       & 
 
  296                             uabs, z1(i,j), pbl(i,j),      & 
 
  297                             z0m(i,j), z0h(i,j), z0e(i,j), & 
 
  298                             dustar, dtstar, dqstar,       & 
 
  299                             dwstar, drlmo, dra,           & 
 
  300                             dfracu10, dfract2, dfracq2    ) 
 
  302              emis = ( 1.0_rp - albedo(i,j,i_r_diffuse,i_r_ir) ) * stb * tmps1(i,j)**4
 
  304              lwd  = rflxd(i,j,i_r_diffuse,i_r_ir)
 
  305              lwu  = rflxd(i,j,i_r_diffuse,i_r_ir)  * albedo(i,j,i_r_diffuse,i_r_ir) + emis
 
  306              swd  = rflxd(i,j,i_r_direct ,i_r_nir) &
 
  307                   + rflxd(i,j,i_r_diffuse,i_r_nir) &
 
  308                   + rflxd(i,j,i_r_direct ,i_r_vis) &
 
  309                   + rflxd(i,j,i_r_diffuse,i_r_vis)
 
  310              swu  = rflxd(i,j,i_r_direct ,i_r_nir) * albedo(i,j,i_r_direct ,i_r_nir) &
 
  311                   + rflxd(i,j,i_r_diffuse,i_r_nir) * albedo(i,j,i_r_diffuse,i_r_nir) &
 
  312                   + rflxd(i,j,i_r_direct ,i_r_vis) * albedo(i,j,i_r_direct ,i_r_vis) &
 
  313                   + rflxd(i,j,i_r_diffuse,i_r_vis) * albedo(i,j,i_r_diffuse,i_r_vis)
 
  316              flx_qv = min( - rhos(i,j) * ustar(i,j) * qstar(i,j) * ra / ( ra+rb(i,j) ), wstr(i,j)/real(dt,rp) )
 
  317              res = swd - swu + lwd - lwu                         &
 
  318                  + cpdry   * rhos(i,j) * ustar(i,j) * tstar(i,j) &
 
  320                  - tc_dz(i,j) * ( tmps1(i,j) - tg(i,j) )
 
  323              dres = -4.0_rp * emis / tmps1(i,j) &
 
  324                   + cpdry   * rhos(i,j) * ( ustar(i,j)*(dtstar-tstar(i,j))/dts0 + tstar(i,j)*(dustar-ustar(i,j))/dts0 )                       &
 
  325                   + lh(i,j) * rhos(i,j) * ( ustar(i,j)*(dqstar-qstar(i,j))/dts0 + qstar(i,j)*(dustar-ustar(i,j))/dts0 ) * ra / ( ra+rb(i,j) ) &
 
  329              if (      abs(res     ) < cpl_phy_sfc_skin_res_min &
 
  330                   .OR. abs(res/dres) < cpl_phy_sfc_skin_err_min ) 
then 
  335              if ( dres < 0.0_rp ) 
then 
  336                 if ( abs(res) > abs(oldres) ) 
then 
  337                    redf = max( tfa*abs(redf), redf_min )
 
  339                    redf = min( tfb*abs(redf), redf_max )
 
  346              dts = - redf * res / dres
 
  347              dts = sign( min( abs(dts), abs(olddts) ), dts )
 
  348              tmps1(i,j) = tmps1(i,j) + dts
 
  356           tmps1(i,j) = min( max( tmps1(i,j),                                                 &
 
  357                                  tmps(i,j) - cpl_phy_sfc_skin_dts_max * real(dt,kind=rp) ), &
 
  358                                  tmps(i,j) + cpl_phy_sfc_skin_dts_max * real(dt,kind=rp) )
 
  360           if ( n > cpl_phy_sfc_skin_itr_max ) 
then 
  362              log_warn(
"CPL_PHY_SFC_skin",*) 
'surface tempearture was not converged. ', trim(model_name)
 
  364              log_info_cont(
'(A,I32)'   ) 
'number of i                        [no unit]  :', i
 
  365              log_info_cont(
'(A,I32)'   ) 
'number of j                        [no unit]  :', j
 
  367              log_info_cont(
'(A,I32)'   ) 
'loop number                        [no unit]  :', n
 
  368              log_info_cont(
'(A,F32.16)') 
'Residual                           [J/m2/s]   :', res
 
  369              log_info_cont(
'(A,F32.16)') 
'delta Residual                     [J/m2/s]   :', dres
 
  371              log_info_cont(
'(A,F32.16)') 
'temperature                        [K]        :', tmpa(i,j)
 
  372              log_info_cont(
'(A,F32.16)') 
'pressure                           [Pa]       :', prsa(i,j)
 
  373              log_info_cont(
'(A,F32.16)') 
'velocity w                         [m/s]      :', wa(i,j)
 
  374              log_info_cont(
'(A,F32.16)') 
'velocity u                         [m/s]      :', ua(i,j)
 
  375              log_info_cont(
'(A,F32.16)') 
'velocity v                         [m/s]      :', va(i,j)
 
  376              log_info_cont(
'(A,F32.16)') 
'absolute velocity                  [m/s]      :', uabs
 
  377              log_info_cont(
'(A,F32.16)') 
'density                            [kg/m3]    :', rhoa(i,j)
 
  378              log_info_cont(
'(A,F32.16)') 
'water vapor mass ratio             [kg/kg]    :', qva(i,j)
 
  379              log_info_cont(
'(A,F32.16)') 
'cell center height                 [m]        :', z1(i,j)
 
  380              log_info_cont(
'(A,F32.16)') 
'atmospheric mixing layer height    [m]        :', pbl(i,j)
 
  381              log_info_cont(
'(A,F32.16)') 
'pressure at the surface            [Pa]       :', prss(i,j)
 
  382              log_info_cont(
'(A,F32.16)') 
'downward radiation (IR, direct )   [J/m2/s]   :', rflxd(i,j,i_r_direct ,i_r_ir )
 
  383              log_info_cont(
'(A,F32.16)') 
'downward radiation (IR, diffuse)   [J/m2/s]   :', rflxd(i,j,i_r_diffuse,i_r_ir )
 
  384              log_info_cont(
'(A,F32.16)') 
'downward radiation (NIR,direct )   [J/m2/s]   :', rflxd(i,j,i_r_direct ,i_r_nir)
 
  385              log_info_cont(
'(A,F32.16)') 
'downward radiation (NIR,diffuse)   [J/m2/s]   :', rflxd(i,j,i_r_diffuse,i_r_nir)
 
  386              log_info_cont(
'(A,F32.16)') 
'downward radiation (VIS,direct )   [J/m2/s]   :', rflxd(i,j,i_r_direct ,i_r_vis)
 
  387              log_info_cont(
'(A,F32.16)') 
'downward radiation (VIS,diffuse)   [J/m2/s]   :', rflxd(i,j,i_r_diffuse,i_r_vis)
 
  389              log_info_cont(
'(A,F32.16)') 
'soil temperature                   [K]        :', tg(i,j)
 
  390              log_info_cont(
'(A,F32.16)') 
'soil water                         [kg/m2]    :', wstr(i,j)
 
  391              log_info_cont(
'(A,F32.16)') 
'surface temperature                [K]        :', tmps(i,j)
 
  392              log_info_cont(
'(A,F32.16)') 
'surface density                    [kg/m3]    :', rhos(i,j)
 
  393              log_info_cont(
'(A,F32.16)') 
'efficiency of evaporation          [1]        :', qvef(i,j)
 
  394              log_info_cont(
'(A,F32.16)') 
'surface albedo (IR, direct )       [1]        :', albedo(i,j,i_r_direct ,i_r_ir )
 
  395              log_info_cont(
'(A,F32.16)') 
'surface albedo (IR, diffuse)       [1]        :', albedo(i,j,i_r_diffuse,i_r_ir )
 
  396              log_info_cont(
'(A,F32.16)') 
'surface albedo (NIR,direct )       [1]        :', albedo(i,j,i_r_direct ,i_r_nir)
 
  397              log_info_cont(
'(A,F32.16)') 
'surface albedo (NIR,diffuse)       [1]        :', albedo(i,j,i_r_diffuse,i_r_nir)
 
  398              log_info_cont(
'(A,F32.16)') 
'surface albedo (VIS,direct )       [1]        :', albedo(i,j,i_r_direct ,i_r_vis)
 
  399              log_info_cont(
'(A,F32.16)') 
'surface albedo (VIS,diffuse)       [1]        :', albedo(i,j,i_r_diffuse,i_r_vis)
 
  400              log_info_cont(
'(A,F32.16)') 
'latent heat                        [J/kg]     :', lh(i,j)
 
  401              log_info_cont(
'(A,F32.16)') 
'stomata registance                 [1/s]      :', rb(i,j)
 
  402              log_info_cont(
'(A,F32.16)') 
'thermal conductivity / depth       [J/m2/s/K] :', tc_dz(i,j)
 
  403              log_info_cont(
'(A,F32.16)') 
'roughness length for momemtum      [m]        :', z0m(i,j)
 
  404              log_info_cont(
'(A,F32.16)') 
'roughness length for heat          [m]        :', z0h(i,j)
 
  405              log_info_cont(
'(A,F32.16)') 
'roughness length for vapor         [m]        :', z0e(i,j)
 
  406              log_info_cont(
'(A,F32.16)') 
'time step                          [s]        :', dt
 
  408              log_info_cont(
'(A,F32.16)') 
'friction velocity                  [m/s]      :', ustar(i,j)
 
  409              log_info_cont(
'(A,F32.16)') 
'friction potential temperature     [K]        :', tstar(i,j)
 
  410              log_info_cont(
'(A,F32.16)') 
'friction water vapor mass ratio    [kg/kg]    :', qstar(i,j)
 
  411              log_info_cont(
'(A,F32.16)') 
'free convection velocity scale     [m/s]      :', wstar(i,j)
 
  412              log_info_cont(
'(A,F32.16)') 
'd(friction velocity)               [m/s]      :', dustar
 
  413              log_info_cont(
'(A,F32.16)') 
'd(friction potential temperature)  [K]        :', dtstar
 
  414              log_info_cont(
'(A,F32.16)') 
'd(friction water vapor mass ratio) [kg/kg]    :', dqstar
 
  415              log_info_cont(
'(A,F32.16)') 
'd(free convection velocity scale)  [m/s]      :', dwstar
 
  416              log_info_cont(
'(A,F32.16)') 
'next surface temperature           [K]        :', tmps1(i,j)
 
  419              if ( .NOT. ( res > -1.0_rp .OR. res < 1.0_rp ) ) 
then  
  420                 log_error(
"CPL_PHY_SFC_skin",*) 
'NaN is detected for surface temperature. ', trim(model_name)
 
  421                 log_error_cont(
'(A,I32)'   ) 
'number of i                        [no unit]  :', i
 
  422                 log_error_cont(
'(A,I32)'   ) 
'number of j                        [no unit]  :', j
 
  423                 log_error_cont(
'(A,I32)'   ) 
'loop number                        [no unit]  :', n
 
  424                 log_error_cont(
'(A,F32.16)') 
'temperature                        [K]        :', tmpa(i,j)
 
  425                 log_error_cont(
'(A,F32.16)') 
'pressure                           [Pa]       :', prsa(i,j)
 
  426                 log_error_cont(
'(A,F32.16)') 
'velocity w                         [m/s]      :', wa(i,j)
 
  427                 log_error_cont(
'(A,F32.16)') 
'velocity u                         [m/s]      :', ua(i,j)
 
  428                 log_error_cont(
'(A,F32.16)') 
'velocity v                         [m/s]      :', va(i,j)
 
  429                 log_error_cont(
'(A,F32.16)') 
'absolute velocity                  [m/s]      :', uabs
 
  430                 log_error_cont(
'(A,F32.16)') 
'density                            [kg/m3]    :', rhoa(i,j)
 
  431                 log_error_cont(
'(A,F32.16)') 
'water vapor mass ratio             [kg/kg]    :', qva(i,j)
 
  432                 log_error_cont(
'(A,F32.16)') 
'cell center height                 [m]        :', z1(i,j)
 
  433                 log_error_cont(
'(A,F32.16)') 
'atmospheric mixing layer height    [m]        :', pbl(i,j)
 
  434                 log_error_cont(
'(A,F32.16)') 
'pressure at the surface            [Pa]       :', prss(i,j)
 
  435                 log_error_cont(
'(A,F32.16)') 
'downward radiation (IR, direct )   [J/m2/s]   :', rflxd(i,j,i_r_direct ,i_r_ir )
 
  436                 log_error_cont(
'(A,F32.16)') 
'downward radiation (IR, diffuse)   [J/m2/s]   :', rflxd(i,j,i_r_diffuse,i_r_ir )
 
  437                 log_error_cont(
'(A,F32.16)') 
'downward radiation (NIR,direct )   [J/m2/s]   :', rflxd(i,j,i_r_direct ,i_r_nir)
 
  438                 log_error_cont(
'(A,F32.16)') 
'downward radiation (NIR,diffuse)   [J/m2/s]   :', rflxd(i,j,i_r_diffuse,i_r_nir)
 
  439                 log_error_cont(
'(A,F32.16)') 
'downward radiation (VIS,direct )   [J/m2/s]   :', rflxd(i,j,i_r_direct ,i_r_vis)
 
  440                 log_error_cont(
'(A,F32.16)') 
'downward radiation (VIS,diffuse)   [J/m2/s]   :', rflxd(i,j,i_r_diffuse,i_r_vis)
 
  441                 log_error_cont(
'(A,F32.16)') 
'soil temperature                   [K]        :', tg(i,j)
 
  442                 log_error_cont(
'(A,F32.16)') 
'soil water                         [kg/m2]    :', wstr(i,j)
 
  443                 log_error_cont(
'(A,F32.16)') 
'surface temperature                [K]        :', tmps(i,j)
 
  444                 log_error_cont(
'(A,F32.16)') 
'surface density                    [kg/m3]    :', rhos(i,j)
 
  445                 log_error_cont(
'(A,F32.16)') 
'efficiency of evaporation          [1]        :', qvef(i,j)
 
  446                 log_error_cont(
'(A,F32.16)') 
'surface albedo (IR, direct )       [1]        :', albedo(i,j,i_r_direct ,i_r_ir )
 
  447                 log_error_cont(
'(A,F32.16)') 
'surface albedo (IR, diffuse)       [1]        :', albedo(i,j,i_r_diffuse,i_r_ir )
 
  448                 log_error_cont(
'(A,F32.16)') 
'surface albedo (NIR,direct )       [1]        :', albedo(i,j,i_r_direct ,i_r_nir)
 
  449                 log_error_cont(
'(A,F32.16)') 
'surface albedo (NIR,diffuse)       [1]        :', albedo(i,j,i_r_diffuse,i_r_nir)
 
  450                 log_error_cont(
'(A,F32.16)') 
'surface albedo (VIS,direct )       [1]        :', albedo(i,j,i_r_direct ,i_r_vis)
 
  451                 log_error_cont(
'(A,F32.16)') 
'surface albedo (VIS,diffuse)       [1]        :', albedo(i,j,i_r_diffuse,i_r_vis)
 
  452                 log_error_cont(
'(A,F32.16)') 
'latent heat                        [J/kg]     :', lh(i,j)
 
  453                 log_error_cont(
'(A,F32.16)') 
'stomata registance                 [1/s]      :', rb(i,j)
 
  454                 log_error_cont(
'(A,F32.16)') 
'thermal conductivity / depth       [J/m2/s/K] :', tc_dz(i,j)
 
  455                 log_error_cont(
'(A,F32.16)') 
'roughness length for momemtum      [m]        :', z0m(i,j)
 
  456                 log_error_cont(
'(A,F32.16)') 
'roughness length for heat          [m]        :', z0h(i,j)
 
  457                 log_error_cont(
'(A,F32.16)') 
'roughness length for vapor         [m]        :', z0e(i,j)
 
  458                 log_error_cont(
'(A,F32.16)') 
'time step                          [s]        :', dt
 
  459                 log_error_cont(
'(A,F32.16)') 
'friction velocity                  [m/s]      :', ustar(i,j)
 
  460                 log_error_cont(
'(A,F32.16)') 
'friction potential temperature     [K]        :', tstar(i,j)
 
  461                 log_error_cont(
'(A,F32.16)') 
'friction water vapor mass ratio    [kg/kg]    :', qstar(i,j)
 
  462                 log_error_cont(
'(A,F32.16)') 
'free convection velocity scale     [m/s]      :', wstar(i,j)
 
  463                 log_error_cont(
'(A,F32.16)') 
'd(friction velocity)               [m/s]      :', dustar
 
  464                 log_error_cont(
'(A,F32.16)') 
'd(friction potential temperature)  [K]        :', dtstar
 
  465                 log_error_cont(
'(A,F32.16)') 
'd(friction water vapor mass ratio) [kg/kg]    :', dqstar
 
  466                 log_error_cont(
'(A,F32.16)') 
'd(free convection velocity scale)  [m/s]      :', dwstar
 
  467                 log_error_cont(
'(A,F32.16)') 
'next surface temperature           [K]        :', tmps1(i,j)
 
  473           tmps(i,j) = tmps1(i,j)
 
  478           call qsat( tmps(i,j), rhos(i,j), qvsat )
 
  481           qvs(i,j) = ( 1.0_rp-qvef(i,j) ) * qva(i,j) &
 
  482                    + (        qvef(i,j) ) * qvsat
 
  484           call bulkflux( tmpa(i,j), tmps(i,j),                  & 
 
  485                          prsa(i,j), prss(i,j),                  & 
 
  486                          qva(i,j), qvs(i,j),                  & 
 
  487                          uabs, z1(i,j), pbl(i,j),               & 
 
  488                          z0m(i,j), z0h(i,j), z0e(i,j),          & 
 
  489                          ustar(i,j), tstar(i,j), qstar(i,j),    & 
 
  490                          wstar(i,j), rlmo(i,j), ra,             & 
 
  491                          fracu10(i,j), fract2(i,j), fracq2(i,j) ) 
 
  493           if ( uabs < eps ) 
then 
  498              mflux = - rhos(i,j) * ustar(i,j)**2
 
  499              zmflx(i,j) = mflux * wa(i,j) / uabs
 
  500              xmflx(i,j) = mflux * ua(i,j) / uabs
 
  501              ymflx(i,j) = mflux * va(i,j) / uabs
 
  503           shflx(i,j) = -rhos(i,j) * ustar(i,j) * tstar(i,j) * cpdry
 
  504           qvflx(i,j) = min( - rhos(i,j) * ustar(i,j) * qstar(i,j) * ra / ( ra+rb(i,j) ), wstr(i,j)/real(dt,rp) )
 
  506           emis = ( 1.0_rp-albedo(i,j,i_r_diffuse,i_r_ir) ) * stb * tmps(i,j)**4
 
  508           lwd  = rflxd(i,j,i_r_diffuse,i_r_ir)
 
  509           lwu  = rflxd(i,j,i_r_diffuse,i_r_ir)  * albedo(i,j,i_r_diffuse,i_r_ir) + emis
 
  510           swd  = rflxd(i,j,i_r_direct ,i_r_nir) &
 
  511                + rflxd(i,j,i_r_diffuse,i_r_nir) &
 
  512                + rflxd(i,j,i_r_direct ,i_r_vis) &
 
  513                + rflxd(i,j,i_r_diffuse,i_r_vis)
 
  514           swu  = rflxd(i,j,i_r_direct ,i_r_nir) * albedo(i,j,i_r_direct ,i_r_nir) &
 
  515                + rflxd(i,j,i_r_diffuse,i_r_nir) * albedo(i,j,i_r_diffuse,i_r_nir) &
 
  516                + rflxd(i,j,i_r_direct ,i_r_vis) * albedo(i,j,i_r_direct ,i_r_vis) &
 
  517                + rflxd(i,j,i_r_diffuse,i_r_vis) * albedo(i,j,i_r_diffuse,i_r_vis)
 
  519           gflx(i,j) = tc_dz(i,j) * ( tmps(i,j) - tg(i,j) )
 
  521           lhflx(i,j) = qvflx(i,j) * lh(i,j)
 
  525           res = swd - swu + lwd - lwu - shflx(i,j) - lhflx(i,j) - gflx(i,j)
 
  528           gflx(i,j) = gflx(i,j) + res
 
  551     call bulkflux_diagnose_surface( ia, is, ie, ja, js, je, &
 
  553                                     tmpa(:,:), qva(:,:),                   & 
 
  554                                     tmps(:,:), qvs(:,:),                   & 
 
  555                                     z1(:,:), z0m(:,:), z0h(:,:), z0e(:,:), & 
 
  556                                     u10(:,:), v10(:,:), t2(:,:), q2(:,:),  & 
 
  557                                     mask = calc_flag(:,:),                 & 
 
  558                                     fracu10 = fracu10(:,:),                & 
 
  559                                     fract2 = fract2(:,:),                  & 
 
  560                                     fracq2 = fracq2(:,:)                   )