!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! GNU General Public License !! !! !! !! This file is part of the Flexible Modeling System (FMS). !! !! !! !! FMS is free software; you can redistribute it and/or modify it !! !! under the terms of the GNU General Public License as published by !! !! the Free Software Foundation, either version 3 of the License, or !! !! (at your option) any later version. !! !! !! !! FMS is distributed in the hope that it will be useful, !! !! but WITHOUT ANY WARRANTY; without even the implied warranty of !! !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !! !! GNU General Public License for more details. !! !! !! !! You should have received a copy of the GNU General Public License !! !! along with FMS. if not, see: http://www.gnu.org/licenses/gpl.txt !! !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! module surface_flux_mod ! ! GFDL ! ! ! ! Driver program for the calculation of fluxes on the exchange grids. ! ! ! ! ! ! ! ============================================================================ use fms_mod, only: FATAL, close_file, mpp_pe, mpp_root_pe, write_version_number use fms_mod, only: file_exist, check_nml_error, open_namelist_file, stdlog use monin_obukhov_mod, only: mo_drag, mo_profile use sat_vapor_pres_mod, only: escomp, descomp use constants_mod, only: cp_air, hlv, stefan, rdgas, rvgas, grav, vonkarm use mpp_mod, only: input_nml_file use mpp_mod, only: mpp_clock_id, mpp_clock_begin, mpp_clock_end, CLOCK_LOOP implicit none private ! ==== public interface ====================================================== public surface_flux ! ==== end of public interface =============================================== ! ! ! For the calculation of fluxes on the exchange grids. ! ! ! For the calculation of fluxes on the exchange grids. ! ! ! ! Air temp lowest atmospheric level. ! ! ! Mixing ratio at lowest atmospheric level (kg/kg). ! ! ! Zonal wind velocity at lowest atmospheric level. ! ! ! Meridional wind velocity at lowest atmospheric level. ! ! ! Pressure lowest atmospheric level. ! ! ! Height lowest atmospheric level. ! ! ! Pressure at the earth's surface ! ! ! Temp at the earth's surface ! ! ! Air temp at the canopy ! ! ! Mixing ratio at earth surface (kg/kg). ! ! ! Zonal wind velocity at earth surface. ! ! ! Meridional wind velocity at earth surface. ! ! ! Momentum roughness length ! ! ! Heat roughness length ! ! ! ! ! Scale factor used to topographic roughness calculation ! ! ! Gustiness factor ! ! ! Sensible heat flux ! ! ! Evaporative water flux ! ! ! Radiative energy flux ! ! ! Zonal momentum flux ! ! ! Meridional momentum flux ! ! ! Momentum exchange coefficient ! ! ! Heat exchange coefficient ! ! ! Moisture exchange coefficient ! ! ! Absolute wind at the lowest atmospheric level ! ! ! Turbulent velocity scale ! ! ! Turbulent buoyant scale ! ! ! Turbulent moisture scale ! ! ! Sensible heat flux temperature sensitivity ! ! ! Moisture flux temperature sensitivity ! ! ! Moisture flux humidity sensitivity ! ! ! Radiative energy flux temperature sensitivity ! ! ! Derivative of sensible heat flux over temp at the lowest atmos level. ! ! ! Derivative of water vapor flux over temp at the lowest atmos level. ! ! ! Derivative of zonal wind stress w.r.t the lowest level zonal ! wind speed of the atmos ! ! ! Derivative of meridional wind stress w.r.t the lowest level meridional ! wind speed of the atmos ! ! ! Time step (it is not used presently) ! ! ! Indicates where land exists (true if exchange cell is on land). ! ! ! Indicates where liquid ocean water exists ! (true if exchange cell is on liquid ocean water). ! ! ! True where the exchange cell is active. ! interface surface_flux ! module procedure surface_flux_0d module procedure surface_flux_1d module procedure surface_flux_2d end interface ! !----------------------------------------------------------------------- character(len=*), parameter :: version = '$Id: surface_flux.F90,v 20.0 2013/12/13 23:27:45 fms Exp $' character(len=*), parameter :: tagname = '$Name: tikal $' logical :: do_init = .true. real, parameter :: d622 = rdgas/rvgas real, parameter :: d378 = 1.-d622 real, parameter :: hlars = hlv/rvgas real, parameter :: gcp = grav/cp_air real, parameter :: kappa = rdgas/cp_air real :: d608 = d378/d622 ! d608 set to zero at initialization if the use of ! virtual temperatures is turned off in namelist ! ---- namelist with default values ------------------------------------------ ! ! ! If q_atm_in (specific humidity) is negative (because of numerical truncation), ! then override with 0. ! ! ! If true, use virtual potential temp to calculate the stability of the surface layer. ! if false, use potential temp. ! ! ! An alternative formulation for gustiness calculation. ! A minimum bound on the wind speed used influx calculations, with the bound ! equal to gust_const ! ! ! The derivative of surface wind stress w.r.t. the zonal wind and ! meridional wind are approximated by the same tendency. ! ! ! An option to provide capability to run the Manabe Climate form of the ! surface flux (coded for legacy purposes). ! ! ! Constant for alternative gustiness calculation. ! ! ! Minimum gustiness used when alt_gustiness = false. ! ! ! Use NCAR climate model turbulent flux calculation described by ! Large and Yeager, NCAR Technical Document, 2004 ! ! ! Use NCAR climate model turbulent flux calculation described by ! Large and Yeager, NCAR Technical Document, 2004, using the original ! GFDL implementation, which contains a bug in the specification of ! the exchange coefficient for the sensible heat. This option is available ! for legacy purposes, and is not recommended for new experiments. ! ! ! Reduce saturation vapor pressures to account for seawater salinity. ! ! logical :: no_neg_q = .false. ! for backwards compatibility logical :: use_virtual_temp = .true. logical :: alt_gustiness = .false. logical :: old_dtaudv = .false. logical :: use_mixing_ratio = .false. real :: gust_const = 1.0 real :: gust_min = 0.0 logical :: ncar_ocean_flux = .false. logical :: ncar_ocean_flux_orig = .false. ! for backwards compatibility logical :: raoult_sat_vap = .false. logical :: do_simple = .false. namelist /surface_flux_nml/ no_neg_q, & use_virtual_temp, & alt_gustiness, & gust_const, & gust_min, & old_dtaudv, & use_mixing_ratio, & ncar_ocean_flux, & ncar_ocean_flux_orig, & raoult_sat_vap, & do_simple integer :: id_clock_escomp = 0, id_clock_mo_drag = 0 contains ! ============================================================================ ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! subroutine surface_flux_1d ( & t_atm, q_atm_in, u_atm, v_atm, p_atm, z_atm, & p_surf, t_surf, t_ca, q_surf, & u_surf, v_surf, & rough_mom, rough_heat, rough_moist, rough_scale, gust, & #ifdef USE_PIK_ML flux_q_scale, & #endif flux_t, flux_q, flux_r, flux_u, flux_v, & cd_m, cd_t, cd_q, & w_atm, u_star, b_star, q_star, & dhdt_surf, dedt_surf, dedq_surf, drdt_surf, & dhdt_atm, dedq_atm, dtaudu_atm, dtaudv_atm, & dt, land, seawater, avail ) ! ! slm Mar 28 2002 -- remove agument drag_q since it is just cd_q*wind ! ============================================================================ ! ---- arguments ----------------------------------------------------------- logical, intent(in), dimension(:) :: land, seawater, avail real, intent(in), dimension(:) :: & t_atm, q_atm_in, u_atm, v_atm, & p_atm, z_atm, t_ca, & p_surf, t_surf, u_surf, v_surf, & rough_mom, rough_heat, rough_moist, rough_scale, gust #ifdef USE_PIK_ML real, intent(in), dimension(:) :: flux_q_scale #endif real, intent(out), dimension(:) :: & flux_t, flux_q, flux_r, flux_u, flux_v, & dhdt_surf, dedt_surf, dedq_surf, drdt_surf, & dhdt_atm, dedq_atm, dtaudu_atm,dtaudv_atm, & w_atm, u_star, b_star, q_star, & cd_m, cd_t, cd_q real, intent(inout), dimension(:) :: q_surf real, intent(in) :: dt ! ---- local constants ----------------------------------------------------- ! temperature increment and its reciprocal value for comp. of derivatives real, parameter:: del_temp=0.1, del_temp_inv=1.0/del_temp ! ---- local vars ---------------------------------------------------------- real, dimension(size(t_atm(:))) :: & thv_atm, th_atm, tv_atm, thv_surf, & e_sat, e_sat1, q_sat, q_sat1, p_ratio, & t_surf0, t_surf1, u_dif, v_dif, & rho_drag, drag_t, drag_m, drag_q, rho, & q_atm, q_surf0, dw_atmdu, dw_atmdv, w_gust integer :: i, nbad if (do_init) call surface_flux_init !---- use local value of surf temp ---- t_surf0 = 200. ! avoids out-of-bounds in es lookup where (avail) where (land) t_surf0 = t_ca elsewhere t_surf0 = t_surf endwhere endwhere t_surf1 = t_surf0 + del_temp call mpp_clock_begin(id_clock_escomp) call escomp ( t_surf0, e_sat ) ! saturation vapor pressure call escomp ( t_surf1, e_sat1 ) ! perturbed vapor pressure call mpp_clock_end(id_clock_escomp) if(use_mixing_ratio) then ! surface mixing ratio at saturation q_sat = d622*e_sat /(p_surf-e_sat ) q_sat1 = d622*e_sat1/(p_surf-e_sat1) elseif(do_simple) then !rif:(09/02/09) q_sat = d622*e_sat / p_surf q_sat1 = d622*e_sat1/ p_surf else ! surface specific humidity at saturation q_sat = d622*e_sat /(p_surf-d378*e_sat ) q_sat1 = d622*e_sat1/(p_surf-d378*e_sat1) endif ! initilaize surface air humidity according to surface type where (land) q_surf0 = q_surf ! land calculates it elsewhere q_surf0 = q_sat ! everything else assumes saturated sfc humidity endwhere if (raoult_sat_vap) where (seawater) q_surf0 = 0.98 * q_surf0 ! check for negative atmospheric humidities where(avail) q_atm = q_atm_in if(no_neg_q) then where(avail .and. q_atm_in < 0.0) q_atm = 0.0 endif ! generate information needed by monin_obukhov where (avail) p_ratio = (p_surf/p_atm)**kappa tv_atm = t_atm * (1.0 + d608*q_atm) ! virtual temperature th_atm = t_atm * p_ratio ! potential T, using p_surf as refernce thv_atm = tv_atm * p_ratio ! virt. potential T, using p_surf as reference thv_surf= t_surf0 * (1.0 + d608*q_surf0 ) ! surface virtual (potential) T ! thv_surf= t_surf0 ! surface virtual (potential) T -- just for testing tun off the q_surf u_dif = u_surf - u_atm ! velocity components relative to surface v_dif = v_surf - v_atm endwhere if(alt_gustiness) then do i = 1, size(avail) if (.not.avail(i)) cycle w_atm(i) = max(sqrt(u_dif(i)**2 + v_dif(i)**2), gust_const) ! derivatives of surface wind w.r.t. atm. wind components if(w_atm(i) > gust_const) then dw_atmdu(i) = u_dif(i)/w_atm(i) dw_atmdv(i) = v_dif(i)/w_atm(i) else dw_atmdu(i) = 0.0 dw_atmdv(i) = 0.0 endif enddo else if (gust_min > 0.0) then where(avail) w_gust = max(gust,gust_min) ! minimum gustiness end where else where(avail) w_gust = gust end where endif where(avail) w_atm = sqrt(u_dif*u_dif + v_dif*v_dif + w_gust*w_gust) ! derivatives of surface wind w.r.t. atm. wind components dw_atmdu = u_dif/w_atm dw_atmdv = v_dif/w_atm endwhere endif ! monin-obukhov similarity theory call mpp_clock_begin(id_clock_mo_drag) call mo_drag (thv_atm, thv_surf, z_atm, & rough_mom, rough_heat, rough_moist, w_atm, & cd_m, cd_t, cd_q, u_star, b_star, avail ) call mpp_clock_end(id_clock_mo_drag) ! override with ocean fluxes from NCAR calculation if (ncar_ocean_flux .or. ncar_ocean_flux_orig) then call ncar_ocean_fluxes (w_atm, th_atm, t_surf0, q_atm, q_surf0, z_atm, & seawater, cd_m, cd_t, cd_q, u_star, b_star ) end if where (avail) ! scale momentum drag coefficient on orographic roughness cd_m = cd_m*(log(z_atm/rough_mom+1)/log(z_atm/rough_scale+1))**2 ! surface layer drag coefficients drag_t = cd_t * w_atm drag_q = cd_q * w_atm drag_m = cd_m * w_atm ! density rho = p_atm / (rdgas * tv_atm) ! sensible heat flux rho_drag = cp_air * drag_t * rho flux_t = rho_drag * (t_surf0 - th_atm) ! flux of sensible heat (W/m**2) dhdt_surf = rho_drag ! d(sensible heat flux)/d(surface temperature) dhdt_atm = -rho_drag*p_ratio ! d(sensible heat flux)/d(atmos temperature) ! evaporation rho_drag = drag_q * rho flux_q = rho_drag * (q_surf0 - q_atm) ! flux of water vapor (Kg/(m**2 s)) #ifdef USE_PIK_ML flux_q = flux_q * flux_q_scale #endif where (land) dedq_surf = rho_drag dedt_surf = 0 elsewhere dedq_surf = 0 dedt_surf = rho_drag * (q_sat1 - q_sat) *del_temp_inv endwhere dedq_atm = -rho_drag ! d(latent heat flux)/d(atmospheric mixing ratio) q_star = flux_q / (u_star * rho) ! moisture scale ! ask Chris and Steve K if we still want to keep this for diagnostics q_surf = q_atm + flux_q / (rho*cd_q*w_atm) ! surface specific humidity ! upward long wave radiation flux_r = stefan*t_surf**4 ! (W/m**2) drdt_surf = 4*stefan*t_surf**3 ! d(upward longwave)/d(surface temperature) ! stresses rho_drag = drag_m * rho flux_u = rho_drag * u_dif ! zonal component of stress (Nt/m**2) flux_v = rho_drag * v_dif ! meridional component of stress elsewhere ! zero-out un-available data in output only fields flux_t = 0.0 flux_q = 0.0 flux_r = 0.0 flux_u = 0.0 flux_v = 0.0 dhdt_surf = 0.0 dedt_surf = 0.0 dedq_surf = 0.0 drdt_surf = 0.0 dhdt_atm = 0.0 dedq_atm = 0.0 u_star = 0.0 b_star = 0.0 q_star = 0.0 q_surf = 0.0 w_atm = 0.0 endwhere ! calculate d(stress component)/d(atmos wind component) dtaudu_atm = 0.0 dtaudv_atm = 0.0 if (old_dtaudv) then where(avail) dtaudv_atm = -rho_drag dtaudu_atm = -rho_drag endwhere else where(avail) dtaudu_atm = -cd_m*rho*(dw_atmdu*u_dif + w_atm) dtaudv_atm = -cd_m*rho*(dw_atmdv*v_dif + w_atm) endwhere endif end subroutine surface_flux_1d ! !####################################################################### subroutine surface_flux_0d ( & t_atm_0, q_atm_0, u_atm_0, v_atm_0, p_atm_0, z_atm_0, & p_surf_0, t_surf_0, t_ca_0, q_surf_0, & u_surf_0, v_surf_0, & rough_mom_0, rough_heat_0, rough_moist_0, rough_scale_0, gust_0, & #ifdef USE_PIK_ML flux_q_scale_0, & #endif flux_t_0, flux_q_0, flux_r_0, flux_u_0, flux_v_0, & cd_m_0, cd_t_0, cd_q_0, & w_atm_0, u_star_0, b_star_0, q_star_0, & dhdt_surf_0, dedt_surf_0, dedq_surf_0, drdt_surf_0, & dhdt_atm_0, dedq_atm_0, dtaudu_atm_0, dtaudv_atm_0, & dt, land_0, seawater_0, avail_0 ) ! ---- arguments ----------------------------------------------------------- logical, intent(in) :: land_0, seawater_0, avail_0 real, intent(in) :: & t_atm_0, q_atm_0, u_atm_0, v_atm_0, & p_atm_0, z_atm_0, t_ca_0, & p_surf_0, t_surf_0, u_surf_0, v_surf_0, & rough_mom_0, rough_heat_0, rough_moist_0, rough_scale_0, gust_0 #ifdef USE_PIK_ML real, intent(in) :: flux_q_scale_0 #endif real, intent(out) :: & flux_t_0, flux_q_0, flux_r_0, flux_u_0, flux_v_0, & dhdt_surf_0, dedt_surf_0, dedq_surf_0, drdt_surf_0, & dhdt_atm_0, dedq_atm_0, dtaudu_atm_0,dtaudv_atm_0, & w_atm_0, u_star_0, b_star_0, q_star_0, & cd_m_0, cd_t_0, cd_q_0 real, intent(inout) :: q_surf_0 real, intent(in) :: dt ! ---- local vars ---------------------------------------------------------- logical, dimension(1) :: land, seawater, avail real, dimension(1) :: & t_atm, q_atm, u_atm, v_atm, & p_atm, z_atm, t_ca, & p_surf, t_surf, u_surf, v_surf, & rough_mom, rough_heat, rough_moist, rough_scale, gust #ifdef USE_PIK_ML real, dimension(1) :: flux_q_scale #endif real, dimension(1) :: & flux_t, flux_q, flux_r, flux_u, flux_v, & dhdt_surf, dedt_surf, dedq_surf, drdt_surf, & dhdt_atm, dedq_atm, dtaudu_atm,dtaudv_atm, & w_atm, u_star, b_star, q_star, & cd_m, cd_t, cd_q real, dimension(1) :: q_surf avail = .true. t_atm(1) = t_atm_0 q_atm(1) = q_atm_0 u_atm(1) = u_atm_0 v_atm(1) = v_atm_0 p_atm(1) = p_atm_0 z_atm(1) = z_atm_0 t_ca(1) = t_ca_0 p_surf(1) = p_surf_0 t_surf(1) = t_surf_0 u_surf(1) = u_surf_0 v_surf(1) = v_surf_0 rough_mom(1) = rough_mom_0 rough_heat(1) = rough_heat_0 rough_moist(1) = rough_moist_0 rough_scale(1) = rough_scale_0 gust(1) = gust_0 #ifdef USE_PIK_ML flux_q_scale(1) = flux_q_scale_0 #endif q_surf(1) = q_surf_0 land(1) = land_0 seawater(1) = seawater_0 avail(1) = avail_0 call surface_flux_1d ( & t_atm, q_atm, u_atm, v_atm, p_atm, z_atm, & p_surf, t_surf, t_ca, q_surf, & u_surf, v_surf, & rough_mom, rough_heat, rough_moist, rough_scale, gust, & #ifdef USE_PIK_ML flux_q_scale, & #endif flux_t, flux_q, flux_r, flux_u, flux_v, & cd_m, cd_t, cd_q, & w_atm, u_star, b_star, q_star, & dhdt_surf, dedt_surf, dedq_surf, drdt_surf, & dhdt_atm, dedq_atm, dtaudu_atm, dtaudv_atm, & dt, land, seawater, avail ) flux_t_0 = flux_t(1) flux_q_0 = flux_q(1) flux_r_0 = flux_r(1) flux_u_0 = flux_u(1) flux_v_0 = flux_v(1) dhdt_surf_0 = dhdt_surf(1) dedt_surf_0 = dedt_surf(1) dedq_surf_0 = dedq_surf(1) drdt_surf_0 = drdt_surf(1) dhdt_atm_0 = dhdt_atm(1) dedq_atm_0 = dedq_atm(1) dtaudu_atm_0 = dtaudu_atm(1) dtaudv_atm_0 = dtaudv_atm(1) w_atm_0 = w_atm(1) u_star_0 = u_star(1) b_star_0 = b_star(1) q_star_0 = q_star(1) q_surf_0 = q_surf(1) cd_m_0 = cd_m(1) cd_t_0 = cd_t(1) cd_q_0 = cd_q(1) end subroutine surface_flux_0d subroutine surface_flux_2d ( & t_atm, q_atm_in, u_atm, v_atm, p_atm, z_atm, & p_surf, t_surf, t_ca, q_surf, & u_surf, v_surf, & rough_mom, rough_heat, rough_moist, rough_scale, gust, & #ifdef USE_PIK_ML flux_q_scale, & #endif flux_t, flux_q, flux_r, flux_u, flux_v, & cd_m, cd_t, cd_q, & w_atm, u_star, b_star, q_star, & dhdt_surf, dedt_surf, dedq_surf, drdt_surf, & dhdt_atm, dedq_atm, dtaudu_atm, dtaudv_atm, & dt, land, seawater, avail ) ! ---- arguments ----------------------------------------------------------- logical, intent(in), dimension(:,:) :: land, seawater, avail real, intent(in), dimension(:,:) :: & t_atm, q_atm_in, u_atm, v_atm, & p_atm, z_atm, t_ca, & p_surf, t_surf, u_surf, v_surf, & rough_mom, rough_heat, rough_moist, rough_scale, gust #ifdef USE_PIK_ML real, intent(in), dimension(:,:) :: flux_q_scale #endif real, intent(out), dimension(:,:) :: & flux_t, flux_q, flux_r, flux_u, flux_v, & dhdt_surf, dedt_surf, dedq_surf, drdt_surf, & dhdt_atm, dedq_atm, dtaudu_atm,dtaudv_atm, & w_atm, u_star, b_star, q_star, & cd_m, cd_t, cd_q real, intent(inout), dimension(:,:) :: q_surf real, intent(in) :: dt ! ---- local vars ----------------------------------------------------------- integer :: j do j = 1, size(t_atm,2) call surface_flux_1d ( & t_atm(:,j), q_atm_in(:,j), u_atm(:,j), v_atm(:,j), p_atm(:,j), z_atm(:,j), & p_surf(:,j), t_surf(:,j), t_ca(:,j), q_surf(:,j), & u_surf(:,j), v_surf(:,j), & rough_mom(:,j), rough_heat(:,j), rough_moist(:,j), rough_scale(:,j), gust(:,j), & #ifdef USE_PIK_ML flux_q_scale(:,j), & #endif flux_t(:,j), flux_q(:,j), flux_r(:,j), flux_u(:,j), flux_v(:,j), & cd_m(:,j), cd_t(:,j), cd_q(:,j), & w_atm(:,j), u_star(:,j), b_star(:,j), q_star(:,j), & dhdt_surf(:,j), dedt_surf(:,j), dedq_surf(:,j), drdt_surf(:,j), & dhdt_atm(:,j), dedq_atm(:,j), dtaudu_atm(:,j), dtaudv_atm(:,j), & dt, land(:,j), seawater(:,j), avail(:,j) ) end do end subroutine surface_flux_2d ! ============================================================================ ! Initialization of the surface flux module--reads the nml. ! subroutine surface_flux_init ! ---- local vars ---------------------------------------------------------- integer :: unit, ierr, io ! read namelist #ifdef INTERNAL_FILE_NML read (input_nml_file, surface_flux_nml, iostat=io) ierr = check_nml_error(io,'surface_flux_nml') #else if ( file_exist('input.nml')) then unit = open_namelist_file () ierr=1; do while (ierr /= 0) read (unit, nml=surface_flux_nml, iostat=io, end=10) ierr = check_nml_error(io,'surface_flux_nml') enddo 10 call close_file (unit) endif #endif ! write version number call write_version_number(version, tagname) unit = stdlog() if ( mpp_pe() == mpp_root_pe() ) write (unit, nml=surface_flux_nml) if(.not. use_virtual_temp) d608 = 0.0 do_init = .false. id_clock_escomp = mpp_clock_id('SFC boundary layer sur_flux_escomp', grain=CLOCK_LOOP) id_clock_mo_drag = mpp_clock_id('SFC boundary layer sur_flux_mo_drag', grain=CLOCK_LOOP) end subroutine surface_flux_init !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! Over-ocean fluxes following Large and Yeager (used in NCAR models) ! ! Original code: GFDL.Climate.Model.Info ! Update Jul2007: GFDL.Climate.Model.Info (ch and ce exchange coeff bugfix) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! subroutine ncar_ocean_fluxes (u_del, t, ts, q, qs, z, avail, & cd, ch, ce, ustar, bstar ) real , intent(in) , dimension(:) :: u_del, t, ts, q, qs, z logical, intent(in) , dimension(:) :: avail real , intent(inout), dimension(:) :: cd, ch, ce, ustar, bstar real :: cd_n10, ce_n10, ch_n10, cd_n10_rt ! neutral 10m drag coefficients real :: cd_rt ! full drag coefficients @ z real :: zeta, x2, x, psi_m, psi_h ! stability parameters real :: u, u10, tv, tstar, qstar, z0, xx, stab integer, parameter :: n_itts = 2 integer i, j if(ncar_ocean_flux_orig) then do i=1,size(u_del(:)) if (avail(i)) then tv = t(i)*(1+0.608*q(i)); u = max(u_del(i), 0.5); ! 0.5 m/s floor on wind (undocumented NCAR) u10 = u; ! first guess 10m wind cd_n10 = (2.7/u10+0.142+0.0764*u10)/1e3; ! L-Y eqn. 6a cd_n10_rt = sqrt(cd_n10); ce_n10 = 34.6 *cd_n10_rt/1e3; ! L-Y eqn. 6b stab = 0.5 + sign(0.5,t(i)-ts(i)) ch_n10 = (18.0*stab+32.7*(1-stab))*cd_n10_rt/1e3; ! L-Y eqn. 6c cd(i) = cd_n10; ! first guess for exchange coeff's at z ch(i) = ch_n10; ce(i) = ce_n10; do j=1,n_itts ! Monin-Obukhov iteration cd_rt = sqrt(cd(i)); ustar(i) = cd_rt*u; ! L-Y eqn. 7a tstar = (ch(i)/cd_rt)*(t(i)-ts(i)); ! L-Y eqn. 7b qstar = (ce(i)/cd_rt)*(q(i)-qs(i)); ! L-Y eqn. 7c bstar(i) = grav*(tstar/tv+qstar/(q(i)+1/0.608)); zeta = vonkarm*bstar(i)*z(i)/(ustar(i)*ustar(i)); ! L-Y eqn. 8a zeta = sign( min(abs(zeta),10.0), zeta ); ! undocumented NCAR x2 = sqrt(abs(1-16*zeta)); ! L-Y eqn. 8b x2 = max(x2, 1.0); ! undocumented NCAR x = sqrt(x2); if (zeta > 0) then psi_m = -5*zeta; ! L-Y eqn. 8c psi_h = -5*zeta; ! L-Y eqn. 8c else psi_m = log((1+2*x+x2)*(1+x2)/8)-2*(atan(x)-atan(1.0)); ! L-Y eqn. 8d psi_h = 2*log((1+x2)/2); ! L-Y eqn. 8e end if u10 = u/(1+cd_n10_rt*(log(z(i)/10)-psi_m)/vonkarm); ! L-Y eqn. 9 cd_n10 = (2.7/u10+0.142+0.0764*u10)/1e3; ! L-Y eqn. 6a again cd_n10_rt = sqrt(cd_n10); ce_n10 = 34.6*cd_n10_rt/1e3; ! L-Y eqn. 6b again stab = 0.5 + sign(0.5,zeta) ch_n10 = (18.0*stab+32.7*(1-stab))*cd_n10_rt/1e3; ! L-Y eqn. 6c again z0 = 10*exp(-vonkarm/cd_n10_rt); ! diagnostic xx = (log(z(i)/10)-psi_m)/vonkarm; cd(i) = cd_n10/(1+cd_n10_rt*xx)**2; ! L-Y 10a xx = (log(z(i)/10)-psi_h)/vonkarm; ch(i) = ch_n10/(1+ch_n10*xx/cd_n10_rt)**2; ! 10b (this code is wrong) ce(i) = ce_n10/(1+ce_n10*xx/cd_n10_rt)**2; ! 10c (this code is wrong) end do end if end do else do i=1,size(u_del(:)) if (avail(i)) then tv = t(i)*(1+0.608*q(i)); u = max(u_del(i), 0.5); ! 0.5 m/s floor on wind (undocumented NCAR) u10 = u; ! first guess 10m wind cd_n10 = (2.7/u10+0.142+0.0764*u10)/1e3; ! L-Y eqn. 6a cd_n10_rt = sqrt(cd_n10); ce_n10 = 34.6 *cd_n10_rt/1e3; ! L-Y eqn. 6b stab = 0.5 + sign(0.5,t(i)-ts(i)) ch_n10 = (18.0*stab+32.7*(1-stab))*cd_n10_rt/1e3; ! L-Y eqn. 6c cd(i) = cd_n10; ! first guess for exchange coeff's at z ch(i) = ch_n10; ce(i) = ce_n10; do j=1,n_itts ! Monin-Obukhov iteration cd_rt = sqrt(cd(i)); ustar(i) = cd_rt*u; ! L-Y eqn. 7a tstar = (ch(i)/cd_rt)*(t(i)-ts(i)); ! L-Y eqn. 7b qstar = (ce(i)/cd_rt)*(q(i)-qs(i)); ! L-Y eqn. 7c bstar(i) = grav*(tstar/tv+qstar/(q(i)+1/0.608)); zeta = vonkarm*bstar(i)*z(i)/(ustar(i)*ustar(i)); ! L-Y eqn. 8a zeta = sign( min(abs(zeta),10.0), zeta ); ! undocumented NCAR x2 = sqrt(abs(1-16*zeta)); ! L-Y eqn. 8b x2 = max(x2, 1.0); ! undocumented NCAR x = sqrt(x2); if (zeta > 0) then psi_m = -5*zeta; ! L-Y eqn. 8c psi_h = -5*zeta; ! L-Y eqn. 8c else psi_m = log((1+2*x+x2)*(1+x2)/8)-2*(atan(x)-atan(1.0)); ! L-Y eqn. 8d psi_h = 2*log((1+x2)/2); ! L-Y eqn. 8e end if u10 = u/(1+cd_n10_rt*(log(z(i)/10)-psi_m)/vonkarm); ! L-Y eqn. 9 cd_n10 = (2.7/u10+0.142+0.0764*u10)/1e3; ! L-Y eqn. 6a again cd_n10_rt = sqrt(cd_n10); ce_n10 = 34.6*cd_n10_rt/1e3; ! L-Y eqn. 6b again stab = 0.5 + sign(0.5,zeta) ch_n10 = (18.0*stab+32.7*(1-stab))*cd_n10_rt/1e3; ! L-Y eqn. 6c again z0 = 10*exp(-vonkarm/cd_n10_rt); ! diagnostic xx = (log(z(i)/10)-psi_m)/vonkarm; cd(i) = cd_n10/(1+cd_n10_rt*xx)**2; ! L-Y 10a xx = (log(z(i)/10)-psi_h)/vonkarm; ch(i) = ch_n10/(1+ch_n10*xx/cd_n10_rt)*sqrt(cd(i)/cd_n10) ! 10b (corrected code) ce(i) = ce_n10/(1+ce_n10*xx/cd_n10_rt)*sqrt(cd(i)/cd_n10) ! 10c (corrected code) end do end if end do endif end subroutine ncar_ocean_fluxes end module surface_flux_mod