#if defined(L08_1A) ! *****************************COPYRIGHT******************************* ! (C) Crown copyright Met Office. All rights reserved. ! For further details please refer to the file COPYRIGHT.txt ! which you should have received as part of this distribution. ! *****************************COPYRIGHT******************************* ! SUBROUTINE SF_EXCH------------------------------------------------ ! ! Purpose: Calculate coefficients of turbulent exchange between ! the surface and the lowest atmospheric layer, and ! "explicit" fluxes between the surface and this layer. ! ! Suitable for Single Column use. ! ! Documentation: UM Documentation Paper No 24, section P243. ! See especially sub-section (ix). ! !--------------------------------------------------------------------- ! Arguments :- SUBROUTINE sf_exch ( & row_length,rows,off_x,off_y,halo_i,halo_j, & land_pts,ntiles,land_index, & tile_index,tile_pts,fland,flandg, & ssi_pts,sea_pts,sice_pts, & ssi_index,sea_index,sice_index,fssi,sea_frac,sice_frac, & nsmax,nsnow,ds,hcons_snow, & bq_1,bt_1,canhc_tile,canopy,catch,dzsoil,flake,gc,hcons, & can_model,catch_snow, lq_mix_bl, & ho2r2_orog,ice_fract,snowdepth,snow_tile,pstar,qw_1,radnet, & radnet_tile,sil_orog,smvcst,tile_frac,timestep, & surf_hgt,emis_tile,emis_soil, & tl_1,ti,ts1, & tsnow, & tstar_tile,tstar_land,tstar_sea,tstar_sice,tstar_ssi,z_land, & l_ctile,seasalinityfactor, & tstar,l_aggregate,l_spec_z0,z0m_scm,z0h_scm,l_dust,l_dust_diag, & vfrac_tile,vshr_land,vshr_ssi,zh,ddmfx, & z0_tile,z1_uv,z1_uv_top,z1_tq,z1_tq_top,land_mask, & su10,sv10,sq1p5,st1p5,sfme,sz0heff,ltimer,formdrag,fd_stab_dep, & orog_drag_param,z0msea, & alpha1,alpha1_sice,ashtf_prime,ashtf_prime_tile,cd,ch, & recip_l_mo_sea,cdr10m, & chr1p5m,chr1p5m_sice,e_sea,fme,fqw_1,fqw_tile,epot_tile, & fqw_ice, & ftl_1,ftl_tile,ftl_ice,fraca,h_blend_orog,h_sea,charnock, & rhostar,resfs,resft,rib,rib_tile, & fb_surf,u_s,q1_sd,t1_sd,z0hssi,z0h_tile,z0h_eff, & z0m_gb,z0mssi,z0m_tile,z0m_eff,rho_aresist,aresist,resist_b, & rho_aresist_tile,aresist_tile,resist_b_tile, & r_b_dust,cd_std_dust,u_s_std_tile, & rho_cd_modv1,rhokh_1,rhokh_1_sice,rhokm_1,rhokm_land,rhokm_ssi, & dtstar_tile,dtstar,rhokh_gb,rhokh_mix,anthrop_heat & ) USE c_z0h_z0m, ONLY : z0h_z0m USE c_vkman USE c_rough USE c_densty USE blend_h USE c_perma USE c_r_cp USE c_g USE c_0_dg_c USE dust_param, ONLY : ndiv USE c_kappai, ONLY : kappai,kappas,dzsea,de USE c_lheat #if defined(UM_RUN) USE blopt8a, ONLY : & surfdivz0t & ,fixed_z0t & ,effective_z0 & ,use_correct_ustar & ,on USE bl_option_mod, ONLY : max_stress_grad #else USE blopt8a, ONLY : & surfdivz0t & ,fixed_z0t & ,effective_z0 & ,use_correct_ustar & ,max_stress_grad & ,on #endif USE nstypes, ONLY : & ! imported scalars with intent(in) npft, urban_canyon, urban_roof USE snow_param, ONLY : cansnowtile & ,rho_snow_const USE jules_mod, ONLY: snowdep_surf USE surf_param, ONLY : h_blend_min & ,ls & ,IP_ScrnDecpl2 & ,z0miz & ,z0sice USE switches, ONLY : & i_modiscopt & ,iseaz0t & ,cor_mo_iter & ,iscrntdiag & ,l_snowdep_surf USE urban_param, ONLY : hgt, hwr, ztm, disp, z0m_mat, z0h_z0m_c, & z0h_z0m_rf USE switches_urban, ONLY : & l_urban2t, l_moruses_rough, l_moruses_storage USE parkind1, ONLY: jprb, jpim USE yomhook, ONLY: lhook, dr_hook IMPLICIT NONE INTEGER & row_length & ! IN Number of X points? ,rows & ! IN Number of Y points? ,off_x & ! Size of small halo in i. ,off_y & ! Size of small halo in j. ,halo_i & ! Size of halo in i direction. ,halo_j & ! Size of halo in j direction. ,land_pts & ! IN No of land points being processed. ,ntiles & ! IN Number of land tiles per land point. ,land_index(land_pts) & ! IN Index of land points. ,tile_index(land_pts,ntiles) & ! IN Index of tile points. ,tile_pts(ntiles) & ! IN Number of tile points. ,nsmax & ! IN Maximum number of snow layers ,nsnow(land_pts,ntiles) & ! IN Number of snow layers ,can_model & ! IN Switch for thermal vegetation canopy. ,ssi_pts & ! IN Number of sea and sea-ice points ,sea_pts & ! IN Number of sea points ,sice_pts & ! IN Number of sea-ice points ,ssi_index(row_length*rows) & ! IN Index of sea and sea-ice points ,sea_index(row_length*rows) & ! IN Index of sea points ,sice_index(row_length*rows) & ! IN Index of sea-ice points ,formdrag & ! IN Switch for orographic form drag ,fd_stab_dep ! IN Switch to implement stability ! dependence of orog form drag REAL & bq_1(row_length,rows) & ! IN A buoyancy parameter for lowest atm ! level ("beta-q twiddle"). ,bt_1(row_length,rows) & ! IN A buoyancy parameter for lowest atm ! level ("beta-T twiddle"). ,canhc_tile(land_pts,ntiles) & ! IN Areal heat capacity of canopy for ! land tiles (J/K/m2). ,canopy(land_pts,ntiles) & ! IN Surface water for land tiles ! (kg/m2). ,catch(land_pts,ntiles) & ! IN Surface capacity (max. surface water) ! of land tiles (kg/m2). ,catch_snow(land_pts,ntiles) & ! IN Snow interception capacity of NLT ! tile (kg/m2). ,ds(land_pts,ntiles,nsmax) & ! IN Snow layer thicknesses (m) ,dzsoil & ! IN Soil or land-ice surface layer ! thickness (m). ,flake(land_pts,ntiles) & ! IN Lake fraction. ,gc(land_pts,ntiles) & ! IN "Stomatal" conductance to evaporation ! for land tiles (m/s). ,hcons(land_pts) & ! IN Soil thermal conductivity including ! effects of water and ice (W/m/K). ,hcons_snow(land_pts,ntiles) & ! IN Snow thermal conductivity (W/m/K) ,ho2r2_orog(land_pts) & ! IN Peak to trough height of unresolved ! orography divided by 2SQRT(2) (m). ,orog_drag_param & ! IN Drag coefficient for orographic ! form drag ,fssi(row_length,rows) & ! IN Fraction of gridbox which is ! sea or sea-ice ,sea_frac(ssi_pts) & ! IN Sea fraction ,sice_frac(ssi_pts) & ! IN Sea-ice fraction ,ice_fract(row_length,rows) & ! IN Fraction of gridbox which is sea-ice. ,fland(land_pts) & ! IN Land fraction on land tiles. ,flandg(row_length,rows) & ! IN Land fraction on all tiles. ,snowdepth(land_pts,ntiles) & ! IN Depth of lying snow (m) ,snow_tile(land_pts,ntiles) & ! IN Lying snow on land tiles (kg/m2). ,pstar(row_length,rows) & ! IN Surface pressure (Pascals). ,qw_1(row_length,rows) & ! IN Total water content of lowest ! atmospheric layer (kg per kg air). ,radnet(row_length,rows) & ! IN Sea-ice net surface radiation (W/m2) ,radnet_tile(land_pts,ntiles) & ! IN Land tile net surface radiation (W/m2) ,anthrop_heat(land_pts,ntiles) & ! IN Anthropogenic Urban heat source (W/m2) ,sil_orog(land_pts) & ! IN Silhouette area of unresolved ! orography per unit horizontal area ,smvcst(land_pts) & ! IN Volumetric saturation point ! - zero at land-ice points. ,tile_frac(land_pts,ntiles) & ! IN Tile fractions. ,timestep & ! IN Timestep in seconds for EPDT calc. ,tl_1(row_length,rows) & ! IN Liquid/frozen water temperature for ! lowest atmospheric layer (K). ,ti(row_length,rows) & ! IN Temperature of sea-ice surface layer ! (K) ,ts1(land_pts) & ! IN Temperature of top soil or land-ice ! layer (K) ,tsnow(land_pts,ntiles,nsmax) & ! IN Snow layer temperatures (K) ,tstar_tile(land_pts,ntiles) & ! IN Tile surface temperatures (K). ,tstar_land(row_length,rows) & ! IN Land mean surface temperature (K). ,tstar_sea(row_length,rows) & ! IN Open sea surface temperature (K). ,tstar_sice(row_length,rows) & ! IN Sea-ice surface temperature (K). ,tstar_ssi(row_length,rows) & ! IN Mean sea surface temperature (K). ,z_land(row_length,rows) & ! IN Land height (m). ,tstar(row_length,rows) & ! IN Gridbox Mean Surface Temperature (K) ,vfrac_tile(land_pts,ntiles) & ! IN Fractional canopy coverage for ! land tiles. ,vshr_land(row_length,rows) & ! IN Magnitude of land sfc-to-lowest-level ! wind shear ,vshr_ssi(row_length,rows) & ! IN Mag. of mean sea sfc-to-lowest-level ! wind shear ,zh(row_length,rows) & ! IN Height above surface of top of ! boundary layer (metres). ,ddmfx(row_length,rows) & ! IN Convective downdraught ! mass-flux at cloud base ,z0_tile(land_pts,ntiles) & ! IN Tile roughness lengths (m). ,z1_uv(row_length,rows) & ! IN Height of lowest uv level (m). ,z1_tq(row_length,rows) & ! IN Height of lowest tq level (m). ! Note, if the grid used is staggered in ! the vertical, Z1_UV and Z1_TQ can be ! different. ,surf_hgt(land_pts,ntiles) & ! IN Height of elevated tile above ! mean gridbox surface (m) ,emis_tile(land_pts,ntiles) & ! IN Emissivity for land tiles ,emis_soil(land_pts) & ! IN Emissivity of underlying soil ,charnock ! Charnock parameter for sea surface REAL, INTENT(IN) :: & z1_uv_top(row_length, rows) ! Height of top of lowest uv-layer REAL, INTENT(IN) :: & z1_tq_top(row_length, rows) ! Height of top of lowest Tq-layer REAL & z0h_scm(row_length,rows) & ! IN Namelist input z0h (if >0) ! (if <=0 use Z0HSEA) ! Used in SCM Configurations ,z0m_scm(row_length,rows)! IN Namelist input z0m (if >0) ! (if <=0 use standard Z0MSEA) ! Used in SCM Configurations LOGICAL & land_mask(row_length,rows) & ! IN .TRUE. for land; .FALSE. elsewhere. ,su10 & ! IN STASH flag for 10-metre W wind. ,sv10 & ! IN STASH flag for 10-metre S wind. ,sq1p5 & ! IN STASH flag for 1.5-metre sp humidity. ,st1p5 & ! IN STASH flag for 1.5-metre temperature. ,sfme & ! IN STASH flag for wind mixing energy flux ,sz0heff & ! IN STASH flag for Z0H_EFF ,lq_mix_bl & ! IN TRUE if mixing ratios used in ! boundary layer code ,ltimer & ! IN Logical for TIMER. ,l_dust & ! IN switch for prognostic mineral dust ,l_dust_diag & ! IN switch for diagnostic mineral dust ! lifting ,l_aggregate & ! IN Logical to set aggregate surface schem ,l_ctile & ! IN switch for coastal tiling ,l_spec_z0 ! IN T if using prescribed ! sea surface roughness lengths REAL, INTENT(IN) :: seasalinityfactor ! Factor allowing for the effect of the salinity of ! sea water on the evaporative flux. ! Modified (INOUT) variables. REAL & z0msea(row_length,rows) ! INOUT Sea-surface roughness length for ! momentum (m). F617. ! Output variables. REAL & alpha1(land_pts,ntiles) & ! OUT Gradients of saturated specific ! humidity with respect to temperature ! between the bottom model layer and ! tile surface ,alpha1_sice(row_length,rows) & ! OUT ALPHA1 for sea-ice. ,ashtf_prime(row_length,rows) & ! OUT Adjusted SEB coefficient for sea-ice ,ashtf_prime_tile(land_pts,ntiles) & ! OUT Adjusted SEB coefficient for land ! points ,cd(row_length,rows) & ! OUT Bulk transfer coefficient for ! momentum. ,cd_ssi(row_length,rows) & ! OUT Bulk transfer coefficient for ! momentum over sea mean. ,ch(row_length,rows) & ! OUT Bulk transfer coefficient for heat ! and/or moisture. ,recip_l_mo_sea(row_length,rows) & ! OUT Reciprocal of the Monin-Obukhov ! length for sea/ice points (m^-1). ,ch_ssi(row_length,rows) & ! OUT Bulk transfer coefficient for heat ! and/or moisture over sea mean. ,cdr10m(1-off_x:row_length+off_x,1-off_y:rows+off_y) & ! OUT Reqd for calculation of 10m wind ! (u & v). ! NBB: This is output on the UV-grid, ! but with the first and last rows set ! to a "missing data indicator". ! Sea-ice leads ignored. ,chr1p5m(land_pts,ntiles) & ! OUT Reqd for calculation of 1.5m temp for ! land tiles. ,chr1p5m_sice(row_length,rows) & ! OUT CHR1P5M for sea and sea-ice ! (leads ignored). ,e_sea(row_length,rows) & ! OUT Evaporation from sea times leads ! fraction (kg/m2/s). Zero over land. ,fme(row_length,rows) & ! OUT Wind mixing energy flux (Watts/sq m). ,fqw_1(row_length,rows) & ! OUT "Explicit" surface flux of QW (i.e. ! evaporation), on P-grid (kg/m2/s). ! for whole grid-box ,fqw_tile(land_pts,ntiles) & ! OUT Local FQW_1 for land tiles. ,fqw_ice(row_length,rows) & ! OUT GBM FQW_1 for sea-ice. ,ftl_1(row_length,rows) & ! OUT "Explicit" surface flux of TL = H/CP. ! (sensible heat / CP). grid-box mean ,ftl_tile(land_pts,ntiles) & ! OUT Local FTL_1 for land tiles. ,ftl_ice(row_length,rows) & ! OUT GBM FTL_1 for sea-ice. ,fraca(land_pts,ntiles) & ! OUT Fraction of surface moisture flux ! with only aerodynamic resistance ! for land tiles. ,h_blend_orog(row_length,rows) & ! OUT Blending height for orographic ! roughness ,h_sea(row_length,rows) & ! OUT Surface sensible heat flux over sea ! times leads fraction (W/m2). ! Zero over land. ,rhostar(row_length,rows) & ! OUT Surface air density ,resfs(land_pts,ntiles) & ! OUT Combined soil, stomatal and ! aerodynamic resistance factor for ! fraction 1-FRACA of land tiles ,resft(land_pts,ntiles) & ! OUT Total resistance factor ! FRACA+(1-FRACA)*RESFS for snow-free ! tiles, 1 for snow and land-ice. ,rib(row_length,rows) & ! OUT Mean bulk Richardson number for ! lowest layer ,rib_tile(land_pts,ntiles) & ! OUT RIB for land tiles. ,fb_surf(row_length,rows) & ! OUT Surface flux buoyancy over ! density (m^2/s^3) ,u_s(row_length,rows) & ! OUT Surface friction velocity (m/s) ,q1_sd(row_length,rows) & ! OUT Standard deviation of turbulent ! fluctuations of surface layer ! specific humidity (kg/kg). ,t1_sd(row_length,rows) & ! OUT Standard deviation of turbulent ! fluctuations of surface layer ! temperature (K). ,z0hssi(row_length,rows) & ! OUT Roughness length for heat and ! moisture over sea/sea ice (m) ,z0h_tile(land_pts,ntiles) & ! OUT Tile roughness lengths for heat ! and moisture ,z0h_eff(row_length,rows) & ! OUT Effective roughness length for ! heat, moisture (m) ,z0m_gb(row_length,rows) & ! OUT Gridbox mean roughness length ! for momentum (m) ,z0mssi(row_length,rows) & ! OUT Roughness length for momentum ! over sea/sea ice (m) ,z0m_tile(land_pts,ntiles) & ! OUT Tile roughness lengths for momentum ,z0m_eff(row_length,rows) & ! OUT Effective roughness length for ! momentum ,rho_aresist(row_length,rows) & ! OUT RHOSTAR*CD_STD*VSHR for SCYCLE ,aresist(row_length,rows) & ! OUT 1/(CD_STD*VSHR) for SCYCLE ,resist_b(row_length,rows) & ! OUT (1/CH-1/CD_STD)/VSHR for SCYCLE ,rho_aresist_tile(land_pts,ntiles) & ! OUT RHOSTAR*CD_STD*VSHR on land tiles ,aresist_tile(land_pts,ntiles) & ! OUT 1/(CD_STD*VSHR) on land tiles ,resist_b_tile(land_pts,ntiles) & ! OUT (1/CH-1/CD_STD)/VSHR on land tiles ,r_b_dust(row_length,rows,ndiv) & ! OUT surf layer res for dust ,cd_std_dust(row_length,rows) & ! OUT Bulk transfer coef. for ! momentum, excluding orographic effects ,u_s_std_tile(land_pts,ntiles) ! OUT Surface layer scaling velocity ! for tiles excluding orographic ! form drag (m/s). ! Surface exchange coefficients;passed to subroutine IMPL_CAL REAL & rho_cd_modv1(row_length,rows) & ! ! OUT rhostar*cD*vshr before horizontal ! ! interpolation output as a diagnostic. ,rhokh_1(land_pts,ntiles) & ! ! OUT Surface exchange coefficient for land ! ! tiles. ,rhokh_1_sice(row_length,rows) & ! ! OUT Surface exchange coefficient for sea ! ! or sea-ice. ,rhokm_1(1-off_x:row_length+off_x,1-off_y:rows+off_y) & ! ! OUT For momentum. NB: This is output on ! ! UV-grid, but with the first and last ! ! rows set to "missing data indicator". ,rhokm_land(1-off_x:row_length+off_x,1-off_y:rows+off_y) & ! ! OUT For land momentum. NB: This is output ! ! on UV-grid, but with the first and ! ! last rows set to "missing data". ,rhokm_ssi(1-off_x:row_length+off_x,1-off_y:rows+off_y) & ! ! OUT For mean sea mom. NB: This is output ! ! on UV-grid, but with the first and ! ! last rows set to "missing data". ,rhokpm(land_pts,ntiles) & ! ! OUT Mixing coefficient for land tiles. ,dtstar_tile(land_pts,ntiles) & ! ! OUT Change in TSTAR over timestep for ! ! land tiles ,dtstar(row_length,rows) & ! ! OUT Change is TSTAR over timestep for ! ! sea-ice ,rhokh_gb(row_length,rows) & ! ! OUT Grid-box surface exchange coefficient ,rhokh_mix(row_length,rows) & ! ! OUT Exchange coeffs for moisture. ,epot_tile(land_pts,ntiles) ! OUT EPOT for land tiles. ! External subprograms called. EXTERNAL qsat_mix,sf_orog,sf_resist,sf_rib, & fcdch,sf_flux,stdev1,sf_orog_gb,sfl_int & ,dustresb EXTERNAL timer ! Define local storage. ! (a) Workspace. REAL & qs1(row_length,rows) & ! Sat. specific humidity ! ! qsat(TL_1,PSTAR) ,z0h_gb(row_length,rows) & ! Gridbox mean roughness length ! ! for heat (m) ,rhokm_ssi_nohalo(row_length,rows) & ! ! like RHOKM_SSI, but with no halo ,lh0 ! Latent heat for snow free surface ! ! =LS for sea-ice, =LC otherwise ! Workspace for sea and sea-ice leads REAL & cd_sea(row_length,rows) & ! Drag coefficient ,ch_sea(row_length,rows) & ! Transfer coefficient for heat and ! ! moisture ,qstar_sea(row_length,rows) & ! Surface saturated sp humidity ,rib_sea(row_length,rows) & ! Bulk Richardson number ,z0h_sea(row_length,rows) & ! Roughness length for heat and ! ! moisture transport ,z0m_sea(row_length,rows) & ! Open sea roughness length for ! ! momentum transport. ,db_sea(row_length,rows) & ! Buoyancy difference for sea points ,v_s_sea(row_length,rows) & ! Surface layer scaling velocity ,alpha1_sea(row_length*rows) & ! ALPHA1 for sea ,ashtf_prime_sea(row_length*rows) & ! ! Adjusted SEB coefficient for sea ! ! points ,hcons_sea(row_length*rows) & ! Heat conductivity into sea ,u_s_std_sea(row_length*rows) & ! Surface friction velocity for sea ! ! (dummy variable for sea) ,v_s_std_sea(row_length*rows) & ! Surface layer scaling velocity ! ! for sea excluding orographic ! ! form drag (m/s). ! ! (dummy variable for sea) ,cd_std_sea(row_length*rows) ! Local drag coefficient for calc ! ! of interpolation coefficient ! ! (dummy variable for sea) ! ! for sea points (m/s). ! Workspace for sea-ice and marginal ice zone REAL & cd_ice(row_length,rows) & ! Drag coefficient ,cd_land(row_length,rows) & ! Bulk transfer coefficient for ! ! momentum over land. ,ch_land(row_length,rows) & ! Bulk transfer coefficient for ! ! het and moisture over land. ,cd_miz(row_length,rows) & ! Drag coefficient ,ch_ice(row_length,rows) & ! Transfer coefficient for heat and ! ! moisture ,ch_miz(row_length,rows) & ! Transfer coefficient for heat and ! ! moisture ,qstar_ice(row_length,rows) & ! Surface saturated sp humidity ,rib_ice(row_length,rows) & ! Bulk Richardson number ,rib_miz(row_length,rows) & ! Bulk Richardson number ,z0_ice(row_length,rows) & ! Roughness length. ,z0_miz(row_length,rows) & ! Roughness length. ,db_ice(row_length,rows) & ! Buoyancy difference for sea ice ,v_s_ice(row_length,rows) & ! Surface layer scaling velocity ! ! for sea ice (m/s). ,v_s_miz(row_length,rows) & ! Surface layer scaling velocity ! ! for marginal sea ice (m/s). ,recip_l_mo_ice(row_length,rows) & ! ! Reciprocal of the Monin-Obukhov ! ! length for sea ice (m^-1). ,recip_l_mo_miz(row_length,rows) & ! ! Reciprocal of the Monin-Obukhov ! ! length for marginal sea ice (m^-1). ,rho_aresist_land(row_length,rows) & ! ! Land mean of rho_aresist_tile ,hcons_sice(row_length*rows) & ! Heat conductiviy into sea-ice ,u_s_std_ice(row_length*rows) & ! Surface friction velocity for sea-i ! ! (dummy variable for sea-ice) ,u_s_std_miz(row_length*rows) & ! Surface friction velocity for ! ! marginal sea-ice ! ! (dummy variable for marginal sea-ic ,v_s_std_ice(row_length*rows) & ! Surface layer scaling velocity ! ! for sea-ice excluding orographic ! ! form drag (m/s). ! ! (dummy variable for sea-ice) ,v_s_std_miz(row_length*rows) & ! Surface layer scaling velocity ! ! for marginal sea-ice excluding ! ! orographic form drag (m/s). ! ! (dummy variable for marginal sea-ic ,cd_std_ice(row_length*rows) & ! Local drag coefficient for calc ! ! of interpolation coefficient ! ! (dummy variable for sea-ice) ,cd_std_miz(row_length*rows) & ! Local drag coefficient for calc ! ! of interpolation coefficient ! ! (dummy variable for marginal sea-ic ,epot_sea(row_length*rows) & ! Potential evaporation from ! sea surface ! (dummy variable for sea surface) ,epot_ice(row_length*rows) ! Potential evaporation from sea-ice ! (dummy variable for sea-ice) REAL & z1_tq_sea(row_length,rows) ! ! Height of lowest model level ! ! relative to sea. REAL, ALLOCATABLE :: & z1_tq_top_sea(:,:) ! Top of lowest Tq layer over sea ! Workspace for land tiles REAL & cd_std(land_pts,ntiles) & ! Local drag coefficient for calc ! ! of interpolation coefficient ,cd_tile(land_pts,ntiles) & ! Drag coefficient ,ch_tile(land_pts,ntiles) & ! Transfer coefficient for heat and ! ! moisture ,chn(land_pts,ntiles) & ! Neutral value of CH. ,ce_d_sh(land_pts) & ! Exposure coefficient * diffusivity ! ! of water vapour * Sherwood number ! ! for canopy snow (m^2/s). ,dq(land_pts) & ! Sp humidity difference between ! ! surface and lowest atmospheric lev ,epdt(land_pts) & ! "Potential" Evaporation * Timestep ,fz0(row_length,rows) & ! Aggregation function for Z0. ,fz0h(row_length,rows) & ! Aggregation function for Z0H. ,pstar_land(land_pts) & ! Surface pressure for land points. ,qstar_tile(land_pts,ntiles) & !Surface saturated sp humidity. ,rhokh_can(land_pts,ntiles) & ! Exchange coefficient for canopy ! ! air to surface ,rhokm_1_tile(land_pts,ntiles) & ! ! Momentum exchange coefficient. ,tsurf(land_pts) & ! Surface layer temp (snow or soil) ( ,dzsurf(land_pts) & ! Surface layer thickness ! ! (snow or soil) (m) ,dzssi(row_length*rows) & ! Surface layer thickness ! ! (sea or sea-ice) (m) ,hcons_surf(land_pts) & ! Thermal conductivity ! ! (snow or soil) (W/m/K) ,wind_profile_factor(land_pts,ntiles) & ! ! For transforming effective surface ! ! transfer coefficients to those ! ! excluding form drag. ,z0m_gb_land(land_pts) & ! GBM momentum land roughness length ,z0h_gb_land(land_pts) & ! GBM land roughness length for heat ,z0m_eff_tile(land_pts,ntiles) & ! ! Effective momentum roughness length ,z0_urban(land_pts,ntiles) & ! ! MORUSES - material r.l. for momentum ,db_tile(land_pts,ntiles) & ! Buoyancy difference for surface ! ! tile ,v_s_tile(land_pts,ntiles) & ! Surface layer scaling velocity ! ! for tiles (m/s). ,v_s_std(land_pts,ntiles) & ! ! Surface layer scaling velocity ! ! for tiles excluding orographic ! ! form drag (m/s). ,u_s_iter_tile(land_pts,ntiles) & ! ! Scaling velocity from middle of ! ! MO scheme - picked up in error by ! ! dust code! ,vshr(row_length,rows) & ! ! Level 1 to surface wind shear ,recip_l_mo_tile(land_pts,ntiles) & ! ! Reciprocal of the Monin-Obukhov ! ! length for tiles (m^-1). ,t_elev(land_pts,ntiles) & ! Temperature at elevated height (k) ,q_elev(land_pts,ntiles) & ! Specific humidity at elevated ! ! height (kg per kg air) ,qs1_elev(land_pts,ntiles) & ! Saturated specific humidity at elev ! ! height (kg per kg air) ,scaling_urban(land_pts,ntiles) ! MORUSES: ground heat flux scaling; ! canyon tile only coupled to soil ! dummy arrays required for sea and se-ice to create universal ! routines for all surfaces REAL & array_zero(row_length*rows) & ! Array of zeros ,array_one(row_length*rows) & ! Array of ones ,array_negone(row_length*rows) ! Array of minus ones INTEGER & array_zero_int(row_length*rows) ! Array of zeros ! (b) Scalars. INTEGER & i,j & ! Loop counter (horizontal field index). ,k & ! Loop counter (tile field index). ,l & ! Loop counter (land point field index). ,n & ! Loop counter (tile index). ,idiv & ! Loop counter (dust division). ,jits ! Counter for iteration for Z0H REAL & tau & ! Magnitude of surface wind stress over sea. ,zetam & ! Temporary in calculation of CHN. ,zetah & ! Temporary in calculation of CHN. ,zeta1 & ! Work space ,z0 & ! yet more workspace ,ustr_l & ! Low-wind estimate of friction velocity ,ustr_n & ! Neutral estimate of friction velocity ,tol_ustr_n & ! Tolerance for USTR_N ,tol_ustr_l & ! Tolerance for USTR_L (see below) ,bl_stress_grad ! Stress gradient across boundary layer INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0 INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1 REAL(KIND=jprb) :: zhook_handle IF (lhook) CALL dr_hook('SF_EXCH',zhook_in,zhook_handle) IF (ltimer) THEN ! DEPENDS ON: timer CALL timer('SFEXCH ',3) END IF array_zero(:)=0.0 array_one(:)=1.0 array_negone(:)=-1.0 array_zero_int(:)=0 ! MORUSES Initialise urban roughness array z0_urban(:,:) = 0.0 DO j=1,rows DO i=1,row_length rhostar(i,j) = pstar(i,j) / ( r*tstar(i,j) ) ! ... surface air density from ideal gas equation END DO END DO !----------------------------------------------------------------------- ! Land surface calculations !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! 0. Initialise FTL_TILE and RIB_TILE on all tiles at all points, ! to allow STASH to process these as diagnostics. !----------------------------------------------------------------------- DO n=1,ntiles DO l=1,land_pts ftl_tile(l,n) = 0.0 rib_tile(l,n) = 0.0 z0m_tile(l,n) = 0.0 u_s_std_tile(l,n)=0.0 END DO END DO DO j=1,rows DO i=1,row_length rhokm_land(i,j) = 0. END DO END DO ! Equivalent snowdepth for surface calculations. snowdep_surf(:,:) = snowdepth(:,:) DO n=1,ntiles IF ( (can_model==4).AND.cansnowtile(n).AND.l_snowdep_surf ) THEN DO l=1,land_pts snowdep_surf(l,n)=snow_tile(l,n)/rho_snow_const END DO END IF END DO !----------------------------------------------------------------------- !! 2. Calculate QSAT values required later. !----------------------------------------------------------------------- ! Calculate temeprature and specific humidity for elevation bands ! DEPENDS ON: elevate CALL elevate( & row_length,rows,land_pts,ntiles,tile_pts,land_index,tile_index, & tl_1,qw_1,qs1,pstar,surf_hgt,t_elev,q_elev, & ltimer & ) ! DEPENDS ON: qsat_mix CALL qsat_mix(qs1,tl_1,pstar,row_length*rows,lq_mix_bl) DO l=1,land_pts j=(land_index(l)-1)/row_length + 1 i = land_index(l) - (j-1)*row_length pstar_land(l) = pstar(i,j) END DO DO n=1,ntiles ! DEPENDS ON: qsat_mix CALL qsat_mix(qstar_tile(:,n),tstar_tile(:,n), & pstar_land,land_pts,lq_mix_bl) ! DEPENDS ON: qsat_mix CALL qsat_mix(qs1_elev(:,n),t_elev(:,n), & pstar_land,land_pts,lq_mix_bl) END DO !----------------------------------------------------------------------- !! 3. Calculation of transfer coefficients and surface layer stability !----------------------------------------------------------------------- !----------------------------------------------------------------------- !! 3.1 Calculate neutral roughness lengths !----------------------------------------------------------------------- ! Land tiles ! Z0_TILE contains the appropriate value for land-ice points, but has to ! be modified for snow-cover on non-land-ice points DO n=1,ntiles DO k=1,tile_pts(n) l = tile_index(k,n) IF ( snow_tile(l,n).gt.0. .AND. & smvcst(l) > EPSILON(smvcst(l)) ) THEN z0 = z0_tile(l,n) - 0.1 * snowdep_surf(l,n) zeta1 = MIN( 5.0e-4 , z0_tile(l,n) ) z0m_tile(l,n) = MAX( zeta1 , z0 ) ELSE z0m_tile(l,n) = z0_tile(l,n) END IF ! MORUSES: z0m_tile (z0_tile previously set in sparm) reset with ztm ! (calculated in init_urban) to remove effects of snow. MORUSES does not ! yet contain a parametrisation for snow although snow should not affect ! the behaviour of bluff bodies rather it affects the material roughness ! of the road & roof ! Note: The test on smvcst is used to alter just the roof tile and not the ! land-ice tiles as the tile is shared at the moment. This will not ! be required when flexible tiles are introduced. IF ( .NOT. l_aggregate & .AND. l_moruses_rough & .AND. smvcst(l) > EPSILON(smvcst) & .AND. ( n == urban_canyon .OR. n == urban_roof ) ) THEN ! Snow could be added here to the material roughness length for ! momentum before passing to urbanz0. Only the road & roof should be ! affected as the walls will be essentially snow-free z0_urban(l,n) = z0m_mat z0m_tile(l,n) = ztm(l) END IF z0h_tile(l,n) = z0h_z0m(n)*z0m_tile(l,n) rib_tile(l,n) = 0. db_tile(l,n) = 0. wind_profile_factor(l,n)=1.0 z0m_eff_tile(l,n)=z0m_tile(l,n) END DO END DO ! MORUSES: urbanz0 calculates r.l. for heat and updates z0h_tile, which is ! previously set above. Two calls are required; one for urban_canyon and one ! for urban_roof. It is done this way to remove the need for do loop over tile ! type and to avoid over-writing any land-ice tiles. If l_aggregate roughness ! lengths are already calculated: z0m in sparm & z0h above. IF ( .NOT. l_aggregate ) THEN IF ( l_moruses_rough ) THEN n = urban_canyon DO k = 1,tile_pts(n) l = tile_index(k,n) j = ( land_index(l) - 1 ) / row_length + 1 i = land_index(l) - ( j - 1 ) * row_length ! DEPENDS ON: urbanz0 CALL urbanz0( & n, z1_uv(i,j), z1_tq(i,j), hgt(l), hwr(l), disp(l), & z0m_mat, z0_urban(l,n), ztm(l), & z0h_tile(l,n) ) ! DEPENDS ON: urbanz0 CALL urbanz0( & urban_roof, z1_uv(i,j), z1_tq(i,j), hgt(l), hwr(l), disp(l), & z0m_mat, z0_urban(l,urban_roof), ztm(l), & z0h_tile(l,urban_roof) ) END DO ELSE IF ( l_urban2t .AND. .NOT. l_moruses_rough ) THEN ! Set values for URBAN-2T n = urban_canyon DO k = 1,tile_pts(n) l = tile_index(k,n) z0h_tile(l,n) = z0h_z0m_c * z0m_tile(l,n) z0h_tile(l,urban_roof) = z0h_z0m_rf * z0m_tile(l,urban_roof) END DO END IF END IF ! Calculate orographic effective parameter for neutral conditions ! if using orographic roughness scheme IF(formdrag == effective_z0) THEN DO n=1,ntiles ! DEPENDS ON: sf_orog CALL sf_orog ( & row_length,rows,land_pts,tile_pts(n), & land_index,tile_index(:,n), & fd_stab_dep,orog_drag_param,ltimer, & ho2r2_orog,rib_tile(:,n),sil_orog,z0m_tile(:,n),z1_uv, & wind_profile_factor(:,n),z0m_eff_tile(:,n) & ) END DO END IF !----------------------------------------------------------------------- ! Calculate RESFT with neutral CH and EPDT=0 for use in calculation ! of Richardson number. RESFT=1 for snow and land-ice. !----------------------------------------------------------------------- DO n=1,ntiles DO k=1,tile_pts(n) l = tile_index(k,n) j=(land_index(l)-1)/row_length + 1 i = land_index(l) - (j-1)*row_length zetam = LOG ( (z1_uv(i,j) + z0m_tile(l,n))/z0m_tile(l,n) ) zetah = LOG ( (z1_tq(i,j) + z0m_tile(l,n))/z0h_tile(l,n) ) chn(l,n) = (vkman/zetah)*(vkman/zetam)* & wind_profile_factor(l,n) dq(l) = qw_1(i,j) - qstar_tile(l,n) epdt(l) = 0.0 END DO ! DEPENDS ON: sf_resist CALL sf_resist ( & row_length,rows,land_pts,tile_pts(n), & land_index,tile_index(:,n), & canopy(:,n),catch(:,n),chn(:,n),dq,epdt,flake(:,n),gc(:,n), & snowdep_surf(:,n),vshr_land,fraca(:,n),resfs(:,n),resft(:,n), & ltimer & ) END DO ! RESFT < 1 for snow on canopy if canopy snow model used ! N.B. chn is calculated in the preceding loop over tiles and ! used in the next loop over npft. This works only if the ! first npft tiles are the vegetated ones. IF (.NOT. l_aggregate .AND. can_model.eq.4) THEN DO n=1,npft IF ( cansnowtile(n) ) THEN DO k=1,tile_pts(n) l = tile_index(k,n) j=(land_index(l)-1)/row_length + 1 i = land_index(l) - (j-1)*row_length IF (snow_tile(l,n) .gt. 0.) THEN gc(l,n) = 0.06*snow_tile(l,n)**0.6*catch_snow(l,n)**0.4 & * 2.06e-5*(tm/tstar_tile(l,n))**1.75 & * (1.79 + 3*SQRT(vshr_land(i,j))) & / (2*rho_ice*5e-4**2) fraca(l,n) = 0. resfs(l,n) = gc(l,n)/(gc(l,n) + chn(l,n)*vshr_land(i,j)) resft(l,n) = resfs(l,n) END IF END DO END IF END DO END IF !----------------------------------------------------------------------- ! 3.2 Calculate bulk Richardson number for the lowest model level. !----------------------------------------------------------------------- ! Land tiles DO n=1,ntiles ! DEPENDS ON: sf_rib CALL sf_rib ( & row_length,rows,land_pts,tile_pts(n), & land_index,tile_index(:,n), & bq_1,bt_1,qstar_tile(:,n),q_elev(:,n),resft(:,n),t_elev(:,n), & tstar_tile(:,n),vshr_land,z0h_tile(:,n),z0m_tile(:,n), & z1_tq,z1_uv, & rib_tile(:,n),db_tile(:,n),ltimer & ) END DO !----------------------------------------------------------------------- !! 3.3 Calculate stability corrected effective roughness length. !! Stability correction only applies to land points. !----------------------------------------------------------------------- IF(formdrag == effective_z0) THEN DO n=1,ntiles ! DEPENDS ON: sf_orog CALL sf_orog ( & row_length,rows,land_pts,tile_pts(n), & land_index,tile_index(:,n), & fd_stab_dep,orog_drag_param,ltimer, & ho2r2_orog,rib_tile(:,n),sil_orog,z0m_tile(:,n),z1_uv, & wind_profile_factor(:,n),z0m_eff_tile(:,n) & ) END DO END IF !----------------------------------------------------------------------- !! 3.4 Calculate CD, CH via routine FCDCH. !----------------------------------------------------------------------- ! Land tiles DO n=1,ntiles ! DEPENDS ON: fcdch CALL fcdch ( & row_length,rows,cor_mo_iter,land_pts,tile_pts(n), & tile_index(:,n),land_index, & db_tile(:,n),vshr_land, & z0m_eff_tile(:,n),z0h_tile(:,n),zh, & z1_uv,z1_uv_top,z1_tq,z1_tq_top,wind_profile_factor(:,n), & ddmfx, & cd_tile(:,n),ch_tile(:,n),cd_std(:,n), & v_s_tile(:,n),v_s_std(:,n),recip_l_mo_tile(:,n), & u_s_iter_tile(:,n), & ltimer & ) END DO IF ( cor_mo_iter == use_correct_ustar ) THEN ! Use correct "standard" ustar DO n=1,ntiles DO k=1,tile_pts(n) l = tile_index(k,n) u_s_std_tile(l,n) = v_s_std(l,n) END DO END DO ELSE ! Use ustar from mid-iteration DO n=1,ntiles DO k=1,tile_pts(n) l = tile_index(k,n) u_s_std_tile(l,n) = u_s_iter_tile(l,n) END DO END DO END IF !----------------------------------------------------------------------- !! 4.1 Recalculate RESFT using "true" CH and EPDT for land tiles !----------------------------------------------------------------------- DO n=1,ntiles DO k=1,tile_pts(n) l = tile_index(k,n) j=(land_index(l)-1)/row_length + 1 i = land_index(l) - (j-1)*row_length dq(l) = qw_1(i,j) - qstar_tile(l,n) epdt(l) = - rhostar(i,j)*ch_tile(l,n)*vshr_land(i,j) & *dq(l)*timestep END DO ! DEPENDS ON: sf_resist CALL sf_resist ( & row_length,rows,land_pts,tile_pts(n), & land_index,tile_index(:,n), & canopy(:,n),catch(:,n),ch_tile(:,n),dq,epdt,flake(:,n), & gc(:,n),snowdep_surf(:,n),vshr_land,fraca(:,n), & resfs(:,n),resft(:,n), & ltimer) END DO IF (.NOT. l_aggregate .AND. can_model.eq.4) THEN DO n=1,npft IF ( cansnowtile(n) ) THEN DO k=1,tile_pts(n) l = tile_index(k,n) j=(land_index(l)-1)/row_length + 1 i = land_index(l) - (j-1)*row_length IF (snow_tile(l,n) .gt. 0.) THEN fraca(l,n) = 0. resfs(l,n) = gc(l,n) / & (gc(l,n) + ch_tile(l,n)*vshr_land(i,j)) resft(l,n) = resfs(l,n) END IF END DO END IF END DO END IF !----------------------------------------------------------------------- ! Calculate gridbox-means of transfer coefficients. !----------------------------------------------------------------------- DO j=1,rows DO i=1,row_length cd_land(i,j) = 0. ch_land(i,j) = 0. END DO END DO ! Land tiles DO n=1,ntiles DO k=1,tile_pts(n) l = tile_index(k,n) j=(land_index(l)-1)/row_length + 1 i = land_index(l) - (j-1)*row_length cd_land(i,j) = cd_land(i,j) + tile_frac(l,n)*cd_tile(l,n) ch_land(i,j) = ch_land(i,j) + tile_frac(l,n)*ch_tile(l,n) END DO END DO !----------------------------------------------------------------------- !! 4.3 Calculate the surface exchange coefficients RHOK(*) and ! resistances for use in Sulphur Cycle ! (Note that CD_STD, CH and VSHR should never = 0) ! RHOSTAR * CD * VSHR stored for diagnostic output before ! horizontal interpolation. !----------------------------------------------------------------------- ! Land tiles DO n=1,ntiles DO k=1,tile_pts(n) l = tile_index(k,n) j=(land_index(l)-1)/row_length + 1 i = land_index(l) - (j-1)*row_length rhokm_1_tile(l,n) = rhostar(i,j)*cd_tile(l,n)*vshr_land(i,j) ! ! P243.124 rhokm_land(i,j) = rhokm_land(i,j) + & tile_frac(l,n)*rhokm_1_tile(l,n) rhokh_1(l,n) = rhostar(i,j)*ch_tile(l,n)*vshr_land(i,j) ! ! P243.125 END DO END DO !----------------------------------------------------------------------- !! Calculate local and gridbox-average surface fluxes of heat and !! moisture. !----------------------------------------------------------------------- DO j=1,rows DO i=1,row_length ftl_1(i,j) = 0. fqw_1(i,j) = 0. END DO END DO ! Land tiles DO n=1,ntiles DO l = 1,land_pts ftl_tile(l,n) = 0. fqw_tile(l,n) = 0. END DO END DO ! Adjust ASHTF for sens. heat flux to ground beneath coniferous canopy rhokh_can(:,:)=0.0 IF (.NOT. l_aggregate .AND. can_model.eq.4) THEN DO n=1,npft IF ( cansnowtile(n) ) THEN DO k=1,tile_pts(n) l = tile_index(k,n) j=(land_index(l)-1)/row_length + 1 i = land_index(l) - (j-1)*row_length rhokh_can(l,n) = rhostar(i,j) * cp / & (43. / (SQRT(cd_tile(l,n))*vshr_land(i,j))) END DO END IF END DO END IF ! Initialise scaling_urban to 1.0 so that it only affects urban tiles when ! MORUSES used with no aggregation. scaling_urban(:,:) = 1.0 IF ( .NOT. l_aggregate .AND. l_moruses_storage ) THEN n = urban_canyon DO l = 1, land_pts IF ( tile_frac(l,n) > 0.0 ) THEN scaling_urban(l,n) = & ( tile_frac(l,n) + tile_frac(l,urban_roof) ) / & tile_frac(l,n) END IF END DO END IF DO n=1,ntiles lh0=lc DO l = 1,land_pts j=(land_index(l)-1)/row_length + 1 i = land_index(l) - (j-1)*row_length tsurf(l) = ts1(l) + t_elev(l,n) - tl_1(i,j) dzsurf(l) = dzsoil hcons_surf(l) = hcons(l) IF ((nsmax > 0).AND.(nsnow(l,n) > 0)) THEN tsurf(l) = tsnow(l,n,1) dzsurf(l) = ds(l,n,1) hcons_surf(l) = hcons_snow(l,n) END IF END DO ! DEPENDS ON: sf_flux CALL sf_flux ( & row_length,rows,land_pts,tile_pts(n),flandg, & land_index,tile_index(:,n), & nsnow(:,n),n,canhc_tile(:,n),dzsurf,hcons_surf, & qs1_elev(:,n),qstar_tile(:,n),q_elev(:,n), & radnet_tile(:,n),resft(:,n),rhokh_1(:,n),smvcst, & snowdepth(:,n),tile_frac(:,n),timestep,t_elev(:,n),tsurf, & tstar_tile(:,n),vfrac_tile(:,n),rhokh_can(:,n),z0h_tile(:,n), & z0m_eff_tile(:,n),z1_tq,lh0,emis_tile(:,n),emis_soil, & 1.0,anthrop_heat(:,n),scaling_urban(:,n), & fqw_1,ftl_1, & alpha1(:,n),ashtf_prime_tile(:,n),fqw_tile(:,n), & epot_tile(:,n),ftl_tile(:,n), & dtstar_tile(:,n),ltimer,0.0 & ) END DO !----------------------------------------------------------------------- ! 4.4 Calculate the standard deviations of layer 1 turbulent ! fluctuations of temperature and humidity using approximate ! formulae from first order closure. !----------------------------------------------------------------------- DO j=1,rows DO i=1,row_length q1_sd(i,j) = 0. t1_sd(i,j) = 0. END DO END DO ! Land tiles DO n=1,ntiles ! DEPENDS ON: stdev1 CALL stdev1 ( & row_length,rows,land_pts,tile_pts(n), & land_index,tile_index(:,n),flandg, & bq_1,bt_1,fqw_tile(:,n),ftl_tile(:,n),rhokm_1_tile(:,n), & rhostar,vshr_land,z0m_tile(:,n),z1_tq,tile_frac(:,n), & q1_sd,t1_sd,ltimer & ) END DO !----------------------------------------------------------------------- !! Call SFL_INT to calculate CDR10M and CHR1P5M - interpolation coeffs !! used to calculate screen temperature, humidity and 10m winds. !----------------------------------------------------------------------- DO j=1,rows DO i=1,row_length cdr10m(i,j) = 0. END DO END DO ! Land tiles IF (su10 .OR. sv10 .OR. sq1p5 .OR. st1p5 .OR. & (IScrnTDiag == IP_ScrnDecpl2) ) THEN DO n=1,ntiles ! DEPENDS ON: sfl_int CALL sfl_int ( & row_length,rows,off_x,off_y,land_pts,tile_pts(n), & tile_index(:,n),land_index,flandg, & vshr_land,cd_std(:,n),cd_tile(:,n),ch_tile(:,n), & tile_frac(:,n), & z0m_eff_tile(:,n),z0m_tile(:,n),z0h_tile(:,n), & recip_l_mo_tile(1,n), & v_s_tile(:,n),v_s_std(:,n), & z1_uv,z1_tq,db_tile(:,n), & su10,sv10,st1p5,sq1p5, & cdr10m,chr1p5m(:,n),ltimer & ) END DO END IF !----------------------------------------------------------------------- ! Sea and sea-ice surface calculations !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! Calculate height of lowest model level relative to sea. !----------------------------------------------------------------------- z1_tq_sea = z1_tq IF (l_ctile) THEN WHERE( (flandg(:,:) > 0.0) .AND. (flandg(:,:) < 1.0) ) z1_tq_sea(:,:) = z1_tq(:,:) + z_land(:,:) END WHERE ENDIF IF (i_modiscopt == on) THEN ALLOCATE(z1_tq_top_sea(row_length,rows)) z1_tq_top_sea = z1_tq_top IF (l_ctile) THEN WHERE( (flandg(:,:) > 0.0) .AND. (flandg(:,:) < 1.0) ) z1_tq_top_sea(:,:) = z1_tq_top(:,:) + z_land(:,:) END WHERE ENDIF ELSE ! Default allocation ALLOCATE(z1_tq_top_sea(row_length,rows)) ENDIF ! DEPENDS ON: qsat_mix CALL qsat_mix(qstar_sea,tstar_sea,pstar,row_length*rows & ,lq_mix_bl) ! DEPENDS ON: qsat_mix CALL qsat_mix(qstar_ice,tstar_sice,pstar,row_length*rows & ,lq_mix_bl) ! Sea, sea-ice leads, sea-ice and marginal ice zone hcons_sea(:) = kappas hcons_sice(:) = kappai DO j=1,rows DO i=1,row_length z0h_sea(i,j) = z0hsea z0_miz(i,j) = z0miz z0_ice(i,j) = z0sice rib_sea(i,j) = 0. rib_ice(i,j) = 0. db_sea(i,j) = 0. db_ice(i,j) = 0. rhokm_ssi(i,j) = 0. rhokm_ssi_nohalo(i,j) = 0. cd_ssi(i,j) = 0. ch_ssi(i,j) = 0. ftl_ice(i,j)=0. fqw_ice(i,j)=0. h_sea(i,j)=0. e_sea(i,j)=0. END DO END DO IF (iseaz0t == surfdivz0t) THEN ! Composite formulation for thermal roughness lengths, ! incoporating the smooth aerodynamic limit for low ! wind speeds and a value based on surface divergence ! theory for higher wind speeds. ! The friction velocity is diagnosed in the surface ! transfer scheme, using z0m from the previous time-step. ! z0[T,q] is also required but depends on u_* in this ! scheme. For practical purposes, it is sufficient to ! infer it from u_* determined from z0m, but in general ! a given value of z0m corresponds to two values of ! u_*, so we need to know whether we are on the low or ! high wind speed side of the minimum value of z0m. ! If we are on the high side, z0[T,q] will be inversely ! proportional to z0m, but on the low side it may follow ! this relationship, or be aerodynamically smooth. In ! the smooth case we iterate u_* from the expression for ! the roughness length and then take the maximum of the ! smooth and high-wind expressions for z0[T,q]. An ! iteration for the low-wind value of u_*, USTR_L is ! carried out. This will converge to the correct limit only ! on the low-wind side of the minimum, and the standard ! criterion that the gradient of a fixed-point iteration ! should be less than 1 in modulus gievs a more precise ! boundary, TOL_USTR_L. For consistency with earlier versions ! of the modset, hard-wired values are retained for the ! operational value of Charnock's parameter. An additional ! check is required, since z0m can be large at low or at ! high wind-speeds. This is less precise and a fixed ! value of 0.07 is used to test USTR_N, which was determined ! by inspection of a graph of z0m against u_*: it is ! unlikely that this will need to be changed unless ! Charnock's constant is greatly altered. IF (charnock == 0.018) THEN tol_ustr_l = 0.055 tol_ustr_n = 0.07 ELSE tol_ustr_l = 0.75*(1.54e-6*g/(2.0*charnock))**0.33333 tol_ustr_n = 0.07 END IF DO j=1, rows DO i=1, row_length ! We need to infer u_* from Z0M. IF (vshr_ssi(i,j) > 0.0) THEN ! Compute u_* using neutral stratification. ! stratification. ustr_n = vkman * vshr_ssi(i,j) / & LOG(1.0 + z1_uv(i,j) / z0msea(i,j) ) ! Compute u_* using low wind approximation. ustr_l = 1.54e-06 / z0msea(i,j) - 1.0e-05 ! Since Z0M could be large for low and high u_*, we use ! USTR_N as an additional check on the regime. IF ( (ustr_n < tol_ustr_n) .AND. & (ustr_l < tol_ustr_l) ) THEN ! Iterate u_* for low winds. DO jits=1, 5 ustr_l=1.54e-06/(z0msea(i,j)-(charnock/g)*ustr_l**2) & -1.0e-05 END DO ! Take the maximum of the smooth and high-wind values. ! A lower limit is imposed on the friction velocity to ! allow for the exceptional case of very low winds: the ! value of 10^-5 is the same as the limit for the momentum ! roughness length. z0h_sea(i,j) = MAX( 2.52e-6/(ustr_l+1.0e-05), & 2.56e-9/z0msea(i,j) ) ELSE ! Take the high-wind value, but limit it to the molecular ! mean free path (we should not hit this limit ! in practice). z0h_sea(i,j) = MAX( 2.56e-9/z0msea(i,j), 7.0e-08 ) END IF END IF END DO END DO ELSE IF (iseaz0t == fixed_z0t) THEN ! Use a fixed thermal roughness length. z0h_sea(1:row_length,1:rows) = z0hsea END IF IF ( l_spec_z0 ) THEN DO j = 1, rows DO i = 1, row_length IF (z0h_scm(i,j) > 0.0) THEN ! Set Z0H from SCM namelist ! (if specified) for sea points z0h_sea(i,j) = z0h_scm(i,j) END IF ! Z0H_SCM END DO ! I END DO ! J END IF !----------------------------------------------------------------------- !! 3.2 Calculate bulk Richardson number for the lowest model level. !----------------------------------------------------------------------- ! Sea, sea-ice and sea-ice leads ! DEPENDS ON: sf_rib CALL sf_rib ( & row_length,rows,ssi_pts,sea_pts, & ssi_index,sea_index, & bq_1,bt_1,qstar_sea,qw_1,array_one,tl_1, & tstar_sea,vshr_ssi,z0h_sea,z0msea, & z1_tq_sea,z1_uv, & rib_sea,db_sea,ltimer & ) ! DEPENDS ON: sf_rib CALL sf_rib ( & row_length,rows,ssi_pts,sice_pts, & ssi_index,sice_index, & bq_1,bt_1,qstar_ice,qw_1,array_one,tl_1, & tstar_sice,vshr_ssi,z0_ice,z0_ice, & z1_tq_sea,z1_uv, & rib_ice,db_ice,ltimer & ) !----------------------------------------------------------------------- !! 3.4 Calculate CD, CH via routine FCDCH. !----------------------------------------------------------------------- ! DEPENDS ON: fcdch CALL fcdch ( & row_length,rows,cor_mo_iter,ssi_pts,sice_pts, & sice_index,ssi_index, & db_ice,vshr_ssi, & z0_ice,z0_ice,zh, & z1_uv,z1_uv_top,z1_tq_sea,z1_tq_top_sea,array_one, & ddmfx, & cd_ice,ch_ice,cd_std_ice, & v_s_ice,v_s_std_ice,recip_l_mo_ice, & u_s_std_ice, & ltimer & ) ! DEPENDS ON: fcdch CALL fcdch ( & row_length,rows,cor_mo_iter,ssi_pts,sice_pts, & sice_index,ssi_index, & db_ice,vshr_ssi, & z0_miz,z0_miz,zh, & z1_uv,z1_uv_top,z1_tq_sea,z1_tq_top_sea,array_one, & ddmfx, & cd_miz,ch_miz,cd_std_miz, & v_s_miz,v_s_std_miz,recip_l_mo_miz, & u_s_std_miz, & ltimer & ) ! DEPENDS ON: fcdch CALL fcdch ( & row_length,rows,cor_mo_iter,ssi_pts,sea_pts, & sea_index,ssi_index, & db_sea,vshr_ssi, & z0msea,z0h_sea,zh, & z1_uv,z1_uv_top,z1_tq_sea,z1_tq_top_sea,array_one, & ddmfx, & cd_sea,ch_sea,cd_std_sea, & v_s_sea,v_s_std_sea,recip_l_mo_sea, & u_s_std_sea, & ltimer & ) ! z1_tq_top_sea is no longer required. DEALLOCATE(z1_tq_top_sea) ! Sea and sea-ice DO j=1,rows DO i=1,row_length IF (flandg(i,j).lt.1.0 ) THEN IF ( ice_fract(i,j) .lt. 0.7 ) THEN cd_ssi(i,j) = ( ice_fract(i,j)*cd_miz(i,j) + & (0.7-ice_fract(i,j))*cd_sea(i,j) ) / 0.7 ! P2430.5 ch_ssi(i,j) = ( ice_fract(i,j)*ch_miz(i,j) + & (0.7-ice_fract(i,j))*ch_sea(i,j) ) / 0.7 ! P2430.4 ELSE cd_ssi(i,j) = ( (1.0-ice_fract(i,j))*cd_miz(i,j) + & (ice_fract(i,j)-0.7)*cd_ice(i,j) ) / 0.3 ! P2430.7 ch_ssi(i,j) = ( (1.0-ice_fract(i,j))*ch_miz(i,j) + & (ice_fract(i,j)-0.7)*ch_ice(i,j) ) / 0.3 ! P2430.7 END IF END IF END DO END DO ! Sea and sea-ice DO j=1,rows DO i=1,row_length IF ( flandg(i,j).lt.1.0 ) THEN rhokm_ssi(i,j) = rhostar(i,j)*cd_ssi(i,j)*vshr_ssi(i,j) ! ! P243.124 rhokm_ssi_nohalo(i,j) = rhokm_ssi(i,j) rhokh_1_sice(i,j) = rhostar(i,j)*ch_ssi(i,j)*vshr_ssi(i,j) ! ! P243.125 ELSE rhokm_ssi( i,j) = 0.0 rhokm_ssi_nohalo(i,j) = 0.0 rhokh_1_sice( i,j) = 0.0 END IF END DO END DO ! Sea and sea-ice lh0=lc dzssi(:) = dzsea ! DEPENDS ON: sf_flux CALL sf_flux ( & row_length,rows,ssi_pts,sea_pts,fssi, & ssi_index,sea_index, & array_zero_int,0,array_zero,dzssi,hcons_sea, & qs1,qstar_sea,qw_1, & radnet,array_one,rhokh_1_sice,array_negone,array_zero, & sea_frac,timestep,tl_1,tstar_sea,tstar_sea, & array_zero,array_zero,z0h_sea,z0msea,z1_tq_sea,lh0, & array_one,array_one, & seasalinityfactor,array_zero,array_one, & fqw_1,ftl_1, & alpha1_sea,ashtf_prime_sea,e_sea,epot_sea,h_sea, & dtstar,ltimer,1.0 & ) lh0=ls dzssi(:)=de ! DEPENDS ON: sf_flux CALL sf_flux ( & row_length,rows,ssi_pts,sice_pts,fssi, & ssi_index,sice_index, & array_zero_int,0,array_zero,dzssi,hcons_sice, & qs1,qstar_ice,qw_1, & radnet,array_one,rhokh_1_sice,array_negone,array_zero, & sice_frac,timestep,tl_1,ti,tstar_sice, & array_zero,array_zero,z0_ice,z0_ice,z1_tq_sea,lh0, & array_one,array_one, & 1.0,array_zero,array_one, & fqw_1,ftl_1, & alpha1_sice,ashtf_prime,fqw_ice,epot_ice,ftl_ice, & dtstar,ltimer,0.0 & ) ! DEPENDS ON: stdev1 CALL stdev1 ( & row_length,rows,ssi_pts,sea_pts, & ssi_index,sea_index,fssi, & bq_1,bt_1,e_sea,h_sea,rhokm_ssi_nohalo, & rhostar,vshr_ssi,z0msea,z1_tq_sea,sea_frac, & q1_sd,t1_sd,ltimer & ) ! DEPENDS ON: stdev1 CALL stdev1 ( & row_length,rows,ssi_pts,sice_pts, & ssi_index,sice_index,fssi, & bq_1,bt_1,fqw_ice,ftl_ice,rhokm_ssi_nohalo, & rhostar,vshr_ssi,z0_ice,z1_tq_sea,sice_frac, & q1_sd,t1_sd,ltimer & ) !----------------------------------------------------------------------- !! 4.6 For sea points, calculate the wind mixing energy flux and the !! sea-surface roughness length on the P-grid, using time-level n !! quantities. !----------------------------------------------------------------------- DO j=1,rows DO i=1,row_length fme(i,j) = 0.0 IF (flandg(i,j).lt.1.0) THEN tau = rhokm_ssi(i,j) * vshr_ssi(i,j) ! P243.130 IF (ice_fract(i,j) .gt. 0.0) & tau = rhostar(i,j) * cd_sea(i,j) & * vshr_ssi(i,j) * vshr_ssi(i,j) IF (sfme) & fme(i,j) = (1.0-ice_fract(i,j)) * tau * SQRT(tau/rhosea) ! ! P243.96 ! Limit Z0MSEA to 0.154m for TAU very small IF( rhostar(i,j) > EPSILON(0.0))THEN z0msea(i,j) = 1.54e-6 / (SQRT(tau/rhostar(i,j)) + 1.0e-5) & + (charnock/g) * (tau / rhostar(i,j)) z0msea(i,j) = MAX ( z0hsea , z0msea(i,j) ) ! ... P243.B6 (Charnock formula) ! TAU/RHOSTAR is "mod VS squared", see eqn P243.131 ELSE z0msea(i,j) = z0hsea END IF END IF END DO END DO IF ( l_spec_z0 ) THEN ! Check for prescribed surface roughness lengths specified in SCM ! NAMELIST. If specified in the &INPROF then they will be used ! instead of Model calculated values DO j=1,rows DO i=1,row_length IF ( z0m_scm(i,j) > 0.0 ) THEN ! Set z0m from SCM namelist for sea points z0msea(i,j) = z0m_scm(i,j) END IF IF ( z0h_scm(i,j) > 0.0 ) THEN ! Set z0h from SCM namelist for sea points z0h_sea(i,j) = z0h_scm(i,j) END IF END DO END DO END IF ! The following lines have been commented out since unless accompanied by ! similar changes for cd_std_sea and v_s_std_sea, which do not feature in ! JULES they will cause an inconsistency in the diagnosis of 10-m winds. ! The question of how winds and temperatures should be diagnosed over ! sea ice is being reconsidered. !--!!----------------------------------------------------------------------- !--!! If sea ice is present then set RECIP_L_MO_SEA to its value over !--!! the ice, a long-standing choice for the screen diagnostics. !--!! Note that RECIP_L_MO_SEA is also used in BDY_EXPL2 to diagnose !--!! shear-dominated boundary layer types. To preserve bit !--!! reproducibility when screen diagnostics are switched on, !--!! this change has been moved outside the if-test on stash logicals !--!!----------------------------------------------------------------------- !--!DO j=1,rows !--! DO i=1,row_length !--! IF (flandg(i,j) < 1.0 .AND. ice_fract(i,j) > 0. ) THEN !--! recip_l_mo_sea(i,j) = recip_l_mo_ice(i,j) !--! END IF !--! END DO !--!END DO ! Sea and sea-ice (leads ignored) IF (su10 .OR. sv10 .OR. sq1p5 .OR. st1p5 .OR. & (IScrnTDiag == IP_ScrnDecpl2) ) THEN ! DEPENDS ON: sfl_int CALL sfl_int ( & row_length,rows,off_x,off_y,ssi_pts,sea_pts, & sea_index,ssi_index,fssi, & vshr_ssi,cd_std_sea,cd_sea,ch_sea, & sea_frac, & z0msea,z0msea,z0h_sea, & recip_l_mo_sea, & v_s_sea,v_s_std_sea, & z1_uv,z1_tq_sea,db_sea, & su10,sv10,st1p5,sq1p5, & cdr10m,chr1p5m_sice,ltimer & ) ! DEPENDS ON: sfl_int CALL sfl_int ( & row_length,rows,off_x,off_y,ssi_pts,sice_pts, & sice_index,ssi_index,fssi, & vshr_ssi,cd_std_ice,cd_ice,ch_ice, & sice_frac, & z0_ice,z0_ice,z0_ice, & recip_l_mo_ice, & v_s_ice,v_s_std_ice, & z1_uv,z1_tq_sea,db_ice, & su10,sv10,st1p5,sq1p5, & cdr10m,chr1p5m_sice,ltimer & ) END IF !----------------------------------------------------------------------- ! GBM diagnstic calculations !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! Calculate effective roughness lengths, orographic blending heights ! and gridbox-average Richardson numbers. !----------------------------------------------------------------------- DO j=1,rows DO i=1,row_length fz0(i,j) = 0. fz0h(i,j) = 0. rib(i,j) = 0. h_blend_orog(i,j) = h_blend_min END DO END DO DO j=1,rows DO i=1,row_length rib(i,j) = rib(i,j) + & (1.0-flandg(i,j))*(1.0-ice_fract(i,j))*rib_sea(i,j) fz0(i,j) = fz0(i,j) + (1.0-flandg(i,j))*(1.0-ice_fract(i,j)) & / (LOG(lb/z0msea(i,j))**2) fz0h(i,j) = fz0h(i,j) + & (1.0-flandg(i,j))*(1.0-ice_fract(i,j)) / & (LOG(lb/z0msea(i,j))*LOG(lb/z0hsea)) rib(i,j) = rib(i,j) + & (1.0-flandg(i,j))*ice_fract(i,j)*rib_ice(i,j) fz0(i,j) = fz0(i,j) + (1.0-flandg(i,j))*ice_fract(i,j) / & (LOG(lb/z0_ice(i,j))**2) fz0h(i,j) = fz0h(i,j) + (1.0-flandg(i,j))*ice_fract(i,j) / & (LOG(lb/z0_ice(i,j))*LOG(lb/z0_ice(i,j))) END DO END DO DO n=1,ntiles DO k=1,tile_pts(n) l = tile_index(k,n) j=(land_index(l)-1)/row_length + 1 i = land_index(l) - (j-1)*row_length rib(i,j) = rib(i,j) + flandg(i,j)*tile_frac(l,n)*rib_tile(l,n) fz0(i,j) = fz0(i,j) & + flandg(i,j)*tile_frac(l,n) / (LOG(lb/z0m_tile(l,n))**2) fz0h(i,j) = fz0h(i,j) + flandg(i,j)*tile_frac(l,n) / & (LOG(lb/z0m_tile(l,n))*LOG(lb/z0h_tile(l,n))) END DO END DO DO j=1,rows DO i=1,row_length z0m_gb(i,j) = lb * EXP( - SQRT(1./fz0(i,j)) ) z0h_gb(i,j) = lb * EXP( - SQRT(fz0(i,j)) / fz0h(i,j)) z0m_eff(i,j) = z0m_gb(i,j) z0h_eff(i,j) = z0h_gb(i,j) END DO END DO DO l = 1,land_pts j=(land_index(l)-1)/row_length + 1 i = land_index(l) - (j-1)*row_length z0m_gb_land(l) = z0m_gb(i,j) z0h_gb_land(l) = z0h_gb(i,j) END DO IF(formdrag == effective_z0) THEN ! DEPENDS ON: sf_orog_gb CALL sf_orog_gb( & row_length,rows,land_pts,land_index, & land_mask,fd_stab_dep,orog_drag_param, & ho2r2_orog,rib,sil_orog,z0m_gb_land,z1_uv, & h_blend_orog,z0m_eff,sz0heff,z0h_gb_land,z0h_eff,ltimer & ) END IF !----------------------------------------------------------------------- ! Calculate gridbox-means of transfer coefficients. !----------------------------------------------------------------------- DO j=1,rows DO i=1,row_length cd(i,j) = flandg(i,j)*cd_land(i,j) + & (1.-flandg(i,j))*cd_ssi(i,j) ch(i,j) = flandg(i,j)*ch_land(i,j) + & (1.-flandg(i,j))*ch_ssi(i,j) END DO END DO !----------------------------------------------------------------------- !! Set grid-box surface exchange coefficients !----------------------------------------------------------------------- DO j=1,rows DO i=1,row_length rhokm_1(i,j)= flandg(i,j) * rhokm_land(i,j) + & (1.0-flandg(i,j)) * rhokm_ssi(i,j) rho_cd_modv1(i,j) = rhokm_1(i,j) ! diagnostic required for VAR END DO END DO DO j=1,rows DO i=1,row_length IF( flandg(i,j) .lt. 1.0) THEN rhokh_gb(i,j) = (1.0 - flandg(i,j))*rhokh_1_sice(i,j) ELSE rhokh_gb(i,j) = 0.0 END IF END DO END DO DO n=1,ntiles DO k=1,tile_pts(n) l = tile_index(k,n) j=(land_index(l)-1)/row_length + 1 i = land_index(l) - (j-1)*row_length rhokh_gb(i,j) = rhokh_gb(i,j) & + flandg(i,j)*tile_frac(l,n)*rhokh_1(l,n) END DO END DO !----------------------------------------------------------------------- !! Calculate scaling parameters required for non-local BL scheme !----------------------------------------------------------------------- DO j=1,rows DO i=1,row_length u_s(i,j) = SQRT( & flandg(i,j)*cd_land(i,j)*vshr_land(i,j)*vshr_land(i,j) & +(1.-flandg(i,j))*cd_ssi(i,j)*vshr_ssi(i,j)*vshr_ssi(i,j) & ) ! !------------------------------------------------------ ! ! Limit the explicitly calculated surface stress, ! ! used to scale the non-local parametrizations, ! ! such that the implied stress gradient across the BL ! ! is less than Max_Stress_Grad. ! !------------------------------------------------------ bl_stress_grad = u_s(i,j)*u_s(i,j) / zh(i,j) IF (bl_stress_grad > max_stress_grad) THEN u_s(i,j) = SQRT(zh(i,j) * max_stress_grad) END IF fb_surf(i,j) = g * ( bt_1(i,j)*ftl_1(i,j) + & bq_1(i,j)*fqw_1(i,j) ) / rhostar(i,j) END DO END DO ! Calculate parameters required for aerosols ! DEPENDS ON: sf_aero CALL sf_aero ( & row_length,rows,land_pts,ntiles,land_index,tile_index,tile_pts, & ltimer,l_dust,l_dust_diag, & flandg,tile_frac,pstar,rhostar,tstar,vshr_land,vshr_ssi, & cd_ssi,ch_ssi,cd_std,ch_tile,rhokh_gb, & rho_aresist,aresist,resist_b,rho_aresist_tile,aresist_tile, & resist_b_tile,r_b_dust,cd_std_dust,rhokh_mix & ) ! set these roughness lengths which otherwise are unspecified z0mssi(:,:) = z0msea( :,:) z0hssi(:,:) = z0h_sea(:,:) IF (ltimer) THEN ! DEPENDS ON: timer CALL timer('SFEXCH ',4) END IF IF (lhook) CALL dr_hook('SF_EXCH',zhook_out,zhook_handle) RETURN END #endif