module atmos_tracer_utilities_mod ! ! William Cooke ! ! ! Bruce Wyman ! ! ! This code provides some utility routines for atmospheric tracers in the FMS framework. ! ! ! This module gives utility routines which can be used to provide ! consistent removal mechanisms for atmospheric tracers. ! ! In particular it provides schemes for wet and dry deposiiton that ! can be easily utilized. ! ! use fms_mod, only : lowercase, & write_version_number, & stdlog, & mpp_pe, & mpp_root_pe, & error_mesg, & NOTE, FATAL use time_manager_mod, only : time_type use diag_manager_mod, only : send_data, & register_diag_field use tracer_manager_mod, only : query_method, & get_tracer_names, & get_number_tracers, & MAX_TRACER_FIELDS use field_manager_mod, only : MODEL_ATMOS, parse use horiz_interp_mod, only : horiz_interp_type, horiz_interp_init, & horiz_interp_new, horiz_interp, horiz_interp_del use monin_obukhov_mod, only : mo_profile use constants_mod, only : GRAV, & ! acceleration due to gravity [m/s2] RDGAS, & ! gas constant for dry air [J/kg/deg] vonkarm, & PI, & DENS_H2O, & ! Water density [kg/m3] WTMH2O, & ! Water molecular weight [g/mole] WTMAIR, & ! Air molecular weight [g/mole] AVOGNO ! Avogadro's number use interpolator_mod, only : interpolator, & obtain_interpolator_time_slices, & unset_interpolator_time_flag, & interpolate_type use astronomy_mod, only : universal_time implicit none private !----------------------------------------------------------------------- !----- interfaces ------- public wet_deposition, & dry_deposition, & dry_deposition_time_vary, & dry_deposition_endts, & interp_emiss, & atmos_tracer_utilities_end, & atmos_tracer_utilities_init, & get_wetdep_param, & get_rh, & get_w10m, & get_cldf, & sjl_fillz !---- version number ----- character(len=128) :: version = '$Id: atmos_tracer_utilities.F90,v 20.0 2013/12/13 23:24:13 fms Exp $' character(len=128) :: tagname = '$Name: tikal $' logical :: module_is_initialized = .FALSE. character(len=7), parameter :: mod_name = 'tracers' integer, parameter :: max_tracers = MAX_TRACER_FIELDS !----------------------------------------------------------------------- !--- identification numbers for diagnostic fields and axes ---- integer :: id_tracer_ddep(max_tracers), id_tracer_dvel(max_tracers), & id_tracer_wdep_ls(max_tracers), id_tracer_wdep_cv(max_tracers), & id_tracer_wdep_lsin(max_tracers), id_tracer_wdep_cvin(max_tracers),& id_tracer_wdep_lsbc(max_tracers), id_tracer_wdep_cvbc(max_tracers),& id_tracer_reevap_ls(max_tracers), id_tracer_reevap_cv(max_tracers),& id_tracer_wdep_ls_3d(max_tracers) integer :: id_tracer_ddep_cmip(max_tracers) integer :: id_w10m, id_delm integer :: id_u_star, id_b_star, id_rough_mom, id_z_pbl, & id_mo_length_inv, id_vds character(len=32), dimension(max_tracers) :: tracer_names = ' ' character(len=32), dimension(max_tracers) :: tracer_units = ' ' character(len=128), dimension(max_tracers) :: tracer_longnames = ' ' character(len=32), dimension(max_tracers) :: tracer_wdep_names = ' ' character(len=32), dimension(max_tracers) :: tracer_wdep_units = ' ' character(len=128), dimension(max_tracers) :: tracer_wdep_longnames = ' ' character(len=32), dimension(max_tracers) :: tracer_ddep_names = ' ' character(len=32), dimension(max_tracers) :: tracer_dvel_names = ' ' character(len=32), dimension(max_tracers) :: tracer_ddep_units = ' ' character(len=32), dimension(max_tracers) :: tracer_dvel_units = ' ' character(len=128), dimension(max_tracers) :: tracer_ddep_longnames = ' ' character(len=128), dimension(max_tracers) :: tracer_dvel_longnames = ' ' real, allocatable :: blon_out(:,:), blat_out(:,:) !----------------parameter values for the diagnostic units-------------- real, parameter :: mw_air = WTMAIR/1000. ! Convert from [g/mole] to [kg/mole] real, parameter :: mw_h2o = WTMH2O/1000. ! Convert from [g/mole] to [kg/mole] real, parameter :: twopi = 2*PI type wetdep_type character (len=500) :: scheme, text_in_scheme, control real :: Henry_constant real :: Henry_variable real :: frac_in_cloud real :: frac_in_cloud_snow real :: alpha_r real :: alpha_s logical :: Lwetdep, Lgas, Laerosol, Lice end type wetdep_type type(wetdep_type), dimension(:), allocatable :: Wetdep type drydep_type character (len=500) :: scheme, name, control real :: land_dry_dep_vel real :: sea_dry_dep_vel real :: ice_dry_dep_vel real :: snow_dry_dep_vel real :: vegn_dry_dep_vel logical :: Ldrydep end type drydep_type type(drydep_type), dimension(:), allocatable :: Drydep contains ! ! ###################################################################### ! ! ! ! This is a routine to create and register the dry and wet deposition ! fields of the tracers. ! ! ! This routine creates diagnostic names for dry and wet deposition fields of the tracers. ! It takes the tracer name and appends "ddep" for the dry deposition field and "wdep" for ! the wet deposition field. This names can then be entered in the diag_table for ! diagnostic output of the tracer dry and wet deposition. The module name associated with ! these fields in "tracers". The units of the deposition fields are assumed to be kg/m2/s. ! ! ! ! The longitude corners for the local domain. ! ! ! The latitude corners for the local domain. ! ! ! The axes relating to the tracer array. ! ! ! Model time. ! subroutine atmos_tracer_utilities_init(lonb, latb, mass_axes, Time) ! Routine to initialize the tracer identification numbers. ! This registers the 2D fields for the wet and dry deposition. real, dimension(:,:), intent(in) :: lonb, latb integer, dimension(3), intent(in) :: mass_axes type(time_type), intent(in) :: Time integer :: ntrace character(len=20) :: units ='' ! integer :: n, logunit character(len=128) :: name logical :: flag ! Make local copies of the local domain dimensions for use ! in interp_emiss. allocate ( blon_out(size(lonb,1),size(lonb,2))) allocate ( blat_out(size(latb,1),size(latb,2))) ! allocate ( data_out(size(lonb(:))-1, size(latb(:))-1)) blon_out = lonb blat_out = latb do n = 1, max_tracers write ( tracer_names(n), 100 ) n write ( tracer_longnames(n), 102 ) n tracer_units(n) = 'none' enddo 100 format ('tr',i3.3) 102 format ('tracer ',i3.3) call get_number_tracers(MODEL_ATMOS, num_tracers= ntrace) if (ntrace > 0) then allocate (Wetdep(ntrace)) allocate (Drydep(ntrace)) endif do n = 1, ntrace !--- set tracer tendency names where tracer names have changed --- call get_tracer_names(MODEL_ATMOS,n,tracer_names(n),tracer_longnames(n),tracer_units(n)) write (name,100) n if (trim(tracer_names(n)) /= name) then tracer_ddep_names(n) = trim(tracer_names(n)) //'_ddep' tracer_dvel_names(n) = trim(tracer_names(n)) //'_dvel' tracer_wdep_names(n) = trim(tracer_names(n)) //'_wdep' endif write (name,102) n if (trim(tracer_longnames(n)) /= name) then tracer_wdep_longnames(n) = & trim(tracer_longnames(n)) // ' wet deposition for tracers' tracer_ddep_longnames(n) = & trim(tracer_longnames(n)) // ' dry deposition for tracers' tracer_dvel_longnames(n) = & trim(tracer_longnames(n)) // ' dry deposition velocity for tracers' endif select case (trim(tracer_units(n))) case ('mmr') units = 'kg/m2/s' case ('kg/kg') units = 'kg/m2/s' case ('vmr') units = 'mole/m2/s' case ('mol/mol') units = 'mole/m2/s' case ('mole/mole') units = 'mole/m2/s' case default units = trim(tracer_units(n))//' kg/(m2 s)' call error_mesg('atmos_tracer_utilities_init',& ' Dry dep units set to '//trim(units)//' in atmos_tracer_utilities for '//trim(tracer_names(n)),& NOTE) end select flag = query_method ('wet_deposition',MODEL_ATMOS,n, & Wetdep(n)%text_in_scheme,Wetdep(n)%control) call get_wetdep_param(Wetdep(n)%text_in_scheme, & Wetdep(n)%control,& Wetdep(n)%scheme, & Wetdep(n)%Henry_constant, & Wetdep(n)%Henry_variable, & Wetdep(n)%frac_in_cloud, & Wetdep(n)%frac_in_cloud_snow, & Wetdep(n)%alpha_r, Wetdep(n)%alpha_s, & Wetdep(n)%Lwetdep, Wetdep(n)%Lgas, & Wetdep(n)%Laerosol, Wetdep(n)%Lice ) Drydep(n)%Ldrydep = query_method ('dry_deposition', MODEL_ATMOS,& n,Drydep(n)%name, Drydep(n)%control) call get_drydep_param(Drydep(n)%name,Drydep(n)%control, & Drydep(n)%scheme,Drydep(n)%land_dry_dep_vel, & Drydep(n)%sea_dry_dep_vel) ! When formulation of dry deposition is resolved perhaps use the following? ! call get_drydep_param & ! (Drydep(n)%name, Drydep(n)%control, Drydep(n)%scheme, & ! Drydep(n)%land_dry_dep_vel, Drydep(n)%sea_dry_dep_vel, & ! Drydep(n)%ice_dry_dep_vel, Drydep(n)%snow_dry_dep_vel, & ! Drydep(n)%vegn_dry_dep_vel ) ! Register the dry deposition of the n tracers id_tracer_ddep(n) = register_diag_field ( mod_name, & trim(tracer_ddep_names(n)), mass_axes(1:2), Time, & trim(tracer_ddep_longnames(n)), & trim(units), missing_value=-999. ) id_tracer_ddep_cmip(n) = register_diag_field ( mod_name, & trim(tracer_ddep_names(n))//'_cmip', mass_axes(1:2), Time, & trim(tracer_ddep_longnames(n)), & 'kg/m2/s', missing_value=-999. ) ! Register the dry deposition of the n tracers id_tracer_dvel(n) = register_diag_field ( mod_name, & trim(tracer_dvel_names(n)), mass_axes(1:2), Time, & trim(tracer_dvel_longnames(n)), & 'm/s', missing_value=-999. ) ! Register the wet deposition of the n tracers by large scale clouds id_tracer_wdep_ls(n) = register_diag_field ( mod_name, & trim(tracer_wdep_names(n))//'_ls', mass_axes(1:2), Time, & trim(tracer_wdep_longnames(n))//' in large scale', & trim(units), missing_value=-999. ) ! Register the wet deposition of the n tracers by large scale clouds 3d id_tracer_wdep_ls_3d(n) = register_diag_field ( mod_name, & trim(tracer_wdep_names(n))//'_ls_3d', mass_axes(1:3), Time, & trim(tracer_wdep_longnames(n))//' in large scale 3D', & trim(units), missing_value=-999. ) ! Register the wet deposition of the n tracers by convective clouds id_tracer_wdep_cv(n) = register_diag_field ( mod_name, & trim(tracer_wdep_names(n))//'_cv', mass_axes(1:2), Time, & trim(tracer_wdep_longnames(n))//' in convective scheme', & trim(units), missing_value=-999. ) ! Register in-cloud rainout by large scale clouds id_tracer_wdep_lsin(n) = register_diag_field ( mod_name, & trim(tracer_wdep_names(n))//'_lsin', mass_axes(1:2), Time, & trim(tracer_wdep_longnames(n))//' in_cloud by lscale precip', & trim(units), missing_value=-999. ) ! Register below-cloud washout by large scale clouds id_tracer_wdep_lsbc(n) = register_diag_field ( mod_name, & trim(tracer_wdep_names(n))//'_lsbc', mass_axes(1:2), Time, & trim(tracer_wdep_longnames(n))//' below_cloud by lscale precip',& trim(units), missing_value=-999. ) ! Register in-cloud re-evaporation by large scale clouds id_tracer_reevap_ls(n) = register_diag_field ( mod_name, & trim(tracer_names(n))//'_reevap_ls', mass_axes(1:3), Time, & trim(tracer_longnames(n))//' re-evap by lscale clouds', & trim(units), missing_value=-999. ) ! Register in-cloud rainout by convective clouds ! Register in-cloud rainout of the n tracers by convective clouds id_tracer_wdep_cvin(n) = register_diag_field ( mod_name, & trim(tracer_wdep_names(n))//'_cvin', mass_axes(1:2), Time, & trim(tracer_wdep_longnames(n))//' in_cloud by conv precip', & trim(units), missing_value=-999. ) ! Register below-cloud washout by convective clouds id_tracer_wdep_cvbc(n) = register_diag_field ( mod_name, & trim(tracer_wdep_names(n))//'_cvbc', mass_axes(1:2), Time, & trim(tracer_wdep_longnames(n))//' below_cloud by conv precip', & trim(units), missing_value=-999. ) ! Register re-evaporation by convective clouds id_tracer_reevap_cv(n) = register_diag_field ( mod_name, & trim(tracer_names(n))//'_reevap_cv', mass_axes(1:3), Time, & trim(tracer_longnames(n))//' re-evap by conv precip', & trim(units), missing_value=-999. ) enddo ! Register scaling factor to calculate wind speed at 10 meters id_delm = register_diag_field ( mod_name, & 'delm', mass_axes(1:2),Time, & 'Scaling factor', 'none', & missing_value=-999. ) ! Register the wind speed at 10 meters id_w10m = register_diag_field ( mod_name, & 'w10m', mass_axes(1:2),Time, & 'Wind speed at 10 meters', 'm/s', & missing_value=-999. ) id_u_star = register_diag_field ( mod_name, & 'u_star_atm', mass_axes(1:2), Time, & 'u star', & 'm/s', missing_value=-999. ) id_b_star = register_diag_field ( mod_name, & 'b_star_atm', mass_axes(1:2), Time, & 'b star', & 'm/s2', missing_value=-999. ) id_rough_mom = register_diag_field ( mod_name, & 'rough_mom_atm', mass_axes(1:2), Time, & 'rough length z0', & 'm', missing_value=-999. ) id_z_pbl = register_diag_field ( mod_name, & 'z_pbl_atm', mass_axes(1:2), Time, & 'z pbl', & 'm', missing_value=-999. ) id_mo_length_inv = register_diag_field ( mod_name, & 'mo_length_inv_atm', mass_axes(1:2), Time, & 'monin Obukhov length', & 'm', missing_value=-999. ) id_vds = register_diag_field ( mod_name, & 'vds_atm', mass_axes(1:2), Time, & 'vds', & 'm/s', missing_value=-999. ) call write_version_number(version, tagname) if ( mpp_pe() == mpp_root_pe() ) then logunit=stdlog() call write_namelist_values (logunit,ntrace) endif module_is_initialized = .TRUE. end subroutine atmos_tracer_utilities_init !#################################################################### subroutine dry_deposition_time_vary (drydep_data, Time) type(time_type), intent(in) :: Time type(interpolate_type), dimension(:), intent(inout) :: drydep_data integer :: n do n=1,size(drydep_data,1) if (Drydep(n)%Ldrydep .and. Drydep(n)%scheme == 'file') then call obtain_interpolator_time_slices (drydep_data(n), Time) endif end do end subroutine dry_deposition_time_vary !#################################################################### subroutine dry_deposition_endts (drydep_data) type(interpolate_type), dimension(:), intent(inout) :: drydep_data integer :: n do n=1, size(drydep_data,1) call unset_interpolator_time_flag (drydep_data(n)) end do end subroutine dry_deposition_endts !#################################################################### ! ! !####################################################################### ! subroutine write_namelist_values (unit, ntrace) integer, intent(in) :: unit, ntrace integer :: n write (unit,10) do n = 1, ntrace write (unit,11) trim(tracer_wdep_names(n)), & trim(tracer_wdep_longnames(n)), & trim(tracer_wdep_units(n)) write (unit,11) trim(tracer_ddep_names(n)), & trim(tracer_ddep_longnames(n)), & trim(tracer_ddep_units(n)) write (unit,11) trim(tracer_dvel_names(n)), & trim(tracer_dvel_longnames(n)), & 'm/s' enddo 10 format (' &TRACER_DIAGNOSTICS_NML', & /,' TRACER: names longnames (units)') 11 format (a16,2x,a,2x,'(',a,')') end subroutine write_namelist_values ! !####################################################################### ! ! subroutine dry_deposition( n, is, js, u, v, T, pwt, pfull, dz, & u_star, landmask, dsinku, tracer, Time, & Time_next, lon, half_day, drydep_data) ! When formulation of dry deposition is resolved perhaps use the following? ! landfr, seaice_cn, snow_area, & ! vegn_cover, vegn_lai, & ! b_star, z_pbl, rough_mom ) ! ! ! Routine to calculate the fraction of tracer to be removed by dry ! deposition. ! ! ! There are three types of dry deposition coded. ! ! 1) Wind driven derived dry deposition velocity. ! ! 2) Fixed dry deposition velocity. ! ! 3) Dry deposition velocities read in from input file ! There are an addition three types of dry deposition coded ! but presently commented out. ! ! 4) Wind driven derived dry deposition velocity, surface dependent. ! ! 5) Wind driven derived dry deposition velocity, surface and boundary ! layer stability dependent. ! ! 6) Fixed dry deposition velocity, difference between land, snow-covered ! land, sea and ice-covered sea. ! ! The theory behind the wind driven dry deposition velocity calculation ! assumes that the deposition can be modeled as a parallel resistance type ! problem. ! ! Total resistance to HNO3-type dry deposition, !

       R = Ra + Rb
!  resisa = aerodynamic resistance
!  resisb = surface resistance (laminar layer + uptake)
!         = 5/u*  [s/cm]        for neutral stability
!      Vd = 1/R
!

! For the fixed dry deposition velocity, there is no change in the ! deposition velocity but the variation of the depth of the surface ! layer implies that there is variation in the amount deposited. ! ! To utilize this section of code add one of the following lines as ! a method for the tracer of interest in the field table. !

! "dry_deposition","wind_driven","surfr=XXX"
!     where XXX is the total resistance defined above.
!
! "dry_deposition","fixed","land=XXX, sea=YYY"
!     where XXX is the dry deposition velocity (m/s) over land
!       and YYY is the dry deposition velocity (m/s) over sea.
!
! "dry_deposition","file","FILENAME.NC"
!     where FILENAME.NC is the NetCDF file name.
!

! ! ! ! ! The tracer number. ! ! ! Start indices for array (computational indices). ! ! ! U wind field. ! ! ! V wind field. ! ! ! Temperature. ! ! ! Pressure differential of half levels. ! ! ! Full pressure levels. ! ! ! Friction velocity. ! ! ! Longitude. ! ! ! Land - sea mask. ! ! ! Dry deposition data interpolated from input file. ! ! ! ! The amount of tracer in the surface layer which is dry deposited per second. ! ! integer, intent(in) :: n, is, js real, intent(in), dimension(:,:) :: u, v, T, pwt, pfull, u_star, tracer, dz real, intent(in), dimension(:,:) :: lon, half_day logical, intent(in), dimension(:,:) :: landmask ! When formulation of dry deposition is resolved perhaps use the following? !real, intent(in), dimension(:,:) :: landfr, z_pbl, b_star, rough_mom !real, intent(in), dimension(:,:) :: seaice_cn, snow_area, vegn_cover, & ! vegn_lai type(time_type), intent(in) :: Time, Time_next type(interpolate_type),intent(inout) :: drydep_data real, intent(out), dimension(:,:) :: dsinku real,dimension(size(u,1),size(u,2)) :: hwindv,frictv,resisa,drydep_vel !real,dimension(size(u,1),size(u,2)) :: mo_length_inv, vds, rs, k1, k2 integer :: i,j, flagsr, id, jd real :: land_dry_dep_vel, sea_dry_dep_vel, ice_dry_dep_vel, & snow_dry_dep_vel, vegn_dry_dep_vel, & surfr, sear, icer, snowr, vegnr real :: diag_scale real :: factor_tmp, gmt, dv_on, dv_off, dayfrac, vd_night, vd_day, loc_angle logical :: used, diurnal integer :: flag_species, flag_diurnal character(len=10) ::units,names character(len=500) :: name,control,scheme, speciesname,dummy ! Default zero dsinku = 0.0 if (.not. Drydep(n)%Ldrydep) return name =Drydep(n)%name control = Drydep(n)%control scheme = Drydep(n)%scheme land_dry_dep_vel = Drydep(n)%land_dry_dep_vel sea_dry_dep_vel = Drydep(n)%sea_dry_dep_vel !ice_dry_dep_vel = Drydep(n)%ice_dry_dep_vel !snow_dry_dep_vel = Drydep(n)%snow_dry_dep_vel !vegn_dry_dep_vel = Drydep(n)%vegn_dry_dep_vel ! delta z = dp/(rho * grav) ! delta z = RT/g*dp/p pwt = dp/g !dz(:,:) = pwt(:,:)*rdgas*T(:,:)/pfull(:,:) id=size(pfull,1); jd=size(pfull,2) select case(lowercase(scheme)) case('wind_driven') ! Calculate horizontal wind velocity and aerodynamic resistance: ! where xxfm=(u*/u) is drag coefficient, Ra=u/(u*^2), ! and u*=sqrt(momentum flux) is friction velocity. ! !**** Compute dry sinks (loss frequency, need modification when !**** different vdep values are to be used for species) flagsr=parse(control,'surfr',surfr) if(flagsr == 0) surfr=500. hwindv=sqrt(u**2+v**2) frictv=u_star resisa=hwindv/(u_star*u_star) where (frictv .lt. 0.1) frictv=0.1 dsinku = (1./(surfr/frictv + resisa))/dz drydep_vel(:,:) = 0. ! case('sfc_dependent_wind_driven') !! Calculate horizontal wind velocity and aerodynamic resistance: !! where xxfm=(u*/u) is drag coefficient, Ra=u/(u*^2), !! and u*=sqrt(momentum flux) is friction velocity. !! !!**** Compute dry sinks (loss frequency, need modification when !!**** different vdep values are to be used for species) ! flagsr=parse(control,'surfr',surfr) ! if(flagsr == 0) surfr=500. ! ! flagsr=parse(control,'sear',sear) ! if(flagsr == 0) sear=surfr ! ! flagsr=parse(control,'icer',icer) ! if(flagsr == 0) icer=surfr ! ! flagsr=parse(control,'snowr',snowr) ! if(flagsr == 0) snowr=surfr ! ! flagsr=parse(control,'vegnr',vegnr) ! if(flagsr == 0) vegnr=surfr ! ! hwindv=sqrt(u**2+v**2) ! frictv=u_star ! resisa=hwindv/(u_star*u_star) ! where (frictv .lt. 0.1) frictv=0.1 ! drydep_vel(:,:) = (1./(surfr/frictv + resisa))* & ! (landfr(:,:) - snow_area(:,:)) + & ! (1./(snowr/frictv + resisa))*snow_area(:,:) +& ! (1./(sear/frictv + resisa))* & ! (1. - landfr(:,:) - seaice_cn(:,:)) + & ! (1./(icer/frictv + resisa))*seaice_cn(:,:) ! dsinku(:,:) = drydep_vel(:,:) / dz(:,:) ! ! case('sfc_BL_dependent_wind_driven') !! Calculate horizontal wind velocity and aerodynamic resistance: !! where xxfm=(u*/u) is drag coefficient, Ra=u/(u*^2), !! and u*=sqrt(momentum flux) is friction velocity. !! !!**** Compute dry sinks (loss frequency, need modification when !!**** different vdep values are to be used for species) ! flagsr=parse(control,'surfr',surfr) ! if(flagsr == 0) surfr=500. ! ! flagsr=parse(control,'sear',sear) ! if(flagsr == 0) sear=surfr ! ! flagsr=parse(control,'icer',icer) ! if(flagsr == 0) icer=surfr ! ! flagsr=parse(control,'snowr',snowr) ! if(flagsr == 0) snowr=surfr ! ! flagsr=parse(control,'vegnr',vegnr) ! if(flagsr == 0) vegnr=surfr ! ! hwindv=sqrt(u**2+v**2) ! frictv=u_star ! resisa=hwindv/(u_star*u_star) ! where (frictv .lt. 0.1) frictv=0.1 ! mo_length_inv = - vonkarm * b_star/(frictv*frictv) ! where (rough_mom > 0.005) ! k1 = 0.001222 * log10(rough_mom) + 0.003906 ! elsewhere ! k1 = 0.001222 * log10(0.005) + 0.003906 ! endwhere ! where(mo_length_inv < 0) ! k2 = 0.0009 * ( - z_pbl * mo_length_inv)**(2.0/3.0) ! elsewhere ! k2 = 0.0 ! endwhere ! vds = frictv * (k1 + k2) ! rs = 1.0 / vds ! drydep_vel(:,:) = (1./(rs + resisa))* & ! (landfr(:,:) - snow_area(:,:)) + & ! (1./(snowr/frictv + resisa))*snow_area(:,:) + & ! (1./(sear/frictv + resisa))* & ! (1. - landfr(:,:) - seaice_cn(:,:)) + & ! (1./(icer/frictv + resisa))*seaice_cn(:,:) ! dsinku(:,:) = drydep_vel(:,:) / dz(:,:) case('fixed') ! For the moment let's try to calculate the delta-z of the bottom ! layer and using a simple dry deposition velocity times the ! timestep, idt, calculate the fraction of the lowest layer which ! deposits. where (landmask(:,:)) ! dry dep value over the land surface drydep_vel(:,:) = land_dry_dep_vel elsewhere ! dry dep value over the sea surface drydep_vel(:,:) = sea_dry_dep_vel endwhere dsinku(:,:) = drydep_vel(:,:) / dz(:,:) ! case('sfc_dependent_fixed') ! drydep_vel(:,:) = land_dry_dep_vel* & ! (landfr(:,:) - snow_area(:,:)) + & ! snow_dry_dep_vel*snow_area(:,:) + & ! sea_dry_dep_vel* & ! (1. - landfr(:,:) - seaice_cn(:,:)) + & ! ice_dry_dep_vel*seaice_cn(:,:) ! dsinku(:,:) = drydep_vel(:,:) / dz(:,:) case('file') flag_species = parse(control,'name',speciesname) if(flag_species>0) then name = trim(speciesname) else call get_tracer_names(MODEL_ATMOS,n,name) endif flag_diurnal = parse(control,'diurnal',dummy) diurnal = (flag_diurnal > 0) call interpolator( drydep_data, Time, drydep_vel, trim(name), is, js ) if (diurnal) then do j = 1,jd do i = 1,id ! half_day is between 0 and pi, so dv_off btwn 0 to pi, dv_on btwn -pi and 0 dv_off = MIN( 1.2*half_day(i,j), PI ) dv_on = -dv_off dayfrac = dv_off/PI ! apply the mean dep vel during polar day or polar night (or nearby) if (dv_off > 0 .and. dv_off < PI ) then vd_night = MIN(0.001, 0.5*drydep_vel(i,j)) vd_day = ( drydep_vel(i,j)-vd_night*(1.-dayfrac) ) / dayfrac gmt = universal_time(Time) loc_angle = gmt + lon(i,j) - PI if (loc_angle >= PI) loc_angle = loc_angle - twopi if (loc_angle < -PI) loc_angle = loc_angle + twopi if( loc_angle >= dv_off .or. loc_angle <= dv_on ) then drydep_vel(i,j) = vd_night else factor_tmp = loc_angle - dv_on factor_tmp = factor_tmp / MAX(2*dv_off,1.e-6) drydep_vel(i,j) = 0.5*PI*sin(factor_tmp*PI)*(vd_day-vd_night) + vd_night end if end if end do end do end if !(diurnal) dsinku(:,:) = drydep_vel(:,:) / dz(:,:) case('default') drydep_vel(:,:) = 0. end select dsinku(:,:) = MAX(dsinku(:,:), 0.0E+00) where(tracer>0) dsinku=dsinku*tracer elsewhere dsinku=0.0 endwhere ! Now save the dry deposition to the diagnostic manager ! delta z = dp/(rho * grav) ! delta z *rho = dp/g ! tracer(kgtracer/kgair) * dz(m)* rho(kgair/m3) = kgtracer/m2 ! so rho drops out of the equation if (id_tracer_ddep(n) > 0 ) then call get_tracer_names(MODEL_ATMOS,n,names,units=units) select case (trim(units)) case ('vmr') diag_scale = mw_air case ('mol/mol') diag_scale = mw_air case ('mole/mole') diag_scale = mw_air case default diag_scale = 1. end select used = send_data ( id_tracer_ddep(n), dsinku*pwt/diag_scale, Time_next, & is_in =is,js_in=js) endif if (id_tracer_ddep_cmip(n) > 0 ) then call get_tracer_names(MODEL_ATMOS,n,names,units=units) select case (trim(names)) case ('so2') diag_scale = mw_air/0.064 case ('so4') diag_scale = mw_air/0.096 case ('dms') diag_scale = mw_air/0.062 case ('nh3') diag_scale = mw_air/0.017 case default diag_scale = 1. end select used = send_data ( id_tracer_ddep_cmip(n), dsinku*pwt/diag_scale,Time_next, & is_in =is,js_in=js) endif if (id_tracer_dvel(n) > 0 ) then used = send_data ( id_tracer_dvel(n), drydep_vel, Time_next, & is_in =is,js_in=js) end if end subroutine dry_deposition ! ! !####################################################################### ! ! ! subroutine wet_deposition( n, T, pfull, phalf, zfull, zhalf, & rain, snow, qdt, cloud, cloud_frac, & f_snow_berg, rain3d, snow3d, & tracer, tracer_dt, Time, cloud_param, & is, js, dt, sum_wdep_out ) ! ! ! Routine to calculate the fraction of tracer removed by wet deposition ! ! ! ! Tracer number ! ! ! start indices for array (computational indices) ! ! ! Temperature ! ! ! Full level pressure field (Pa) ! ! ! Half level pressure field (Pa) ! ! ! Full level height field (m) ! ! ! Half level height field (m) ! ! ! Precipitation in the form of rain ! ! ! Precipitation in the form of snow ! ! ! The tendency of the specific humidity (+ condenstate) due to the cloud parametrization (kg/kg/s) ! ! ! Cloud amount (liquid + ice) (kg/kg) ! ! ! Cloud area fraction ! ! ! Precipitation in the form of rain (kg/m2/s) ! ! ! Precipitation in the form of snow (kg/m2/s) ! ! ! The tracer field ! ! ! The time structure for submitting wet deposition as a diagnostic ! ! ! Is this a convective (convect) or large scale (lscale) cloud parametrization? ! ! ! The model timestep (in seconds) ! ! ! The tendency of the tracer field due to wet deposition ! ! ! Schemes allowed here are ! ! 1) Deposition removed in the same fractional amount as the modeled precipitation rate is to ! a standardized precipitation rate. ! Basically this scheme assumes that a fractional area of the gridbox is affected by ! precipitation and that this precipitation rate is due to a cloud of standardized cloud ! liquid water content. Removal is constant throughout the column where precipitation is occuring. ! ! 2) Removal according to Henry's Law. This law states that the ratio of the concentation in ! cloud water and the partial pressure in the interstitial air is a constant. If tracer ! is in VMR, the units for Henry's constant are mole/L/Pa (normally it is mole/L/atm). ! Parameters for a large number of species can be found at ! http://www.mpch-mainz.mpg.de/~sander/res/henry.html ! ! 3) Aerosol removal, using specified in-cloud tracer fraction ! 4) Similar as 3) with some lwh modifications ! ! To utilize this section of code add one of the following lines as ! a method for the tracer of interest in the field table. !

! "wet_deposition","henry","henry=XXX, dependence=YYY"
!     where XXX is the Henry's constant for the tracer in question
!       and YYY is the temperature dependence of the Henry's Law constant.
!
! "wet_deposition","fraction","lslwc=XXX, convlwc=YYY"
!     where XXX is the liquid water content of a standard large scale cloud
!       and YYY is the liquid water content of a standard convective cloud.
!

! !----------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------- integer, intent(in) :: n, is, js real, intent(in), dimension(:,:,:) :: T, pfull,phalf, zfull, zhalf, qdt, cloud, tracer real, intent(in), dimension(:,:,:) :: cloud_frac real, intent(in), dimension(:,:,:) :: f_snow_berg ! snow production by Bergeron process real, intent(in), dimension(:,:) :: rain, snow character(len=*), intent(in) :: cloud_param type (time_type), intent(in) :: Time real, intent(out), dimension(:,:,:) :: tracer_dt real, intent(in) :: dt real, intent(in), dimension(:,:,:) :: rain3d, snow3d real, intent(out), dimension(:,:), optional :: sum_wdep_out !----------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------- real, dimension(size(T,1),size(T,2),size(pfull,3)) :: & Htemp, xliq, n_air, rho_air, pwt, zdel, precip3d, scav_fact3d, & precip3ds, precip3dr real, dimension(size(T,1),size(T,2)) :: & temp_factor, scav_factor, washout, sum_wdep, & w_h2o, K1, K2, beta, f_a, scav_factor_s, & wdep_in, wdep_bc, fluxr,fluxs, tracer_flux real, dimension(size(T,1),size(T,2),size(pfull,3)) :: & in_temp, bc_temp, dt_temp, reevap_fraction, reevap_diag integer, dimension(size(T,1),size(T,2)) :: & ktopcd, kendcd real, dimension(size(rain3d,1),size(rain3d,2),size(rain3d,3)) :: rainsnow3d integer :: i, j, k, kk, id, jd, kd, flaglw real, dimension(size(T,3)) :: conc real :: conc_rain, conc_rain_total, conc_sat real, parameter :: DENS_SNOW = 500. ! Snow density [kg/m3] real, parameter :: RGAS = 8.3143 ! ideal gas constant Pa m3/mol/K real :: & Henry_constant, Henry_variable, & clwc, wash, premin, prenow, hwtop, & diag_scale real, parameter :: & inv298p15 = 1./298.15, & ! 1/K kboltz = 1.38E-23, & ! J/K rain_diam = 1.89e-3, & ! mean diameter of rain drop (m) rain_vterm = 7.48, & ! rain drop terminal velocity (m/s) vk_air = 6.18e-6, & ! kinematic viscosity of air (m^2/s) d_g = 1.12e-5, & ! diffusive coefficient (m^2/s) geo_fac = 6., & ! geometry factor (surface area/volume = geo_fac/diameter) cm3_2_m3 = 1.e-6 ! m3/cm3 real :: & k_g, & ! mass transfer coefficient (m/s) stay, & ! fraction fall_time ! fall time through layer (s) real :: f_a0, scav_factor0, sa_drop0, fgas0 real :: frac_in_cloud, frac_in_cloud_snow, frac_int, ph real , parameter :: & R_r = 0.001, & ! radius of cloud-droplets for rain R_s = 0.001, & ! radius of cloud-droplets for snow frac_int_gas = 1.0, & frac_int_aerosol= 0.5 real :: alpha_r, alpha_s logical :: & used, & Lwetdep, Lgas, Laerosol, Lice character(len=500) :: & tracer_name, control, scheme, units, & text_in_scheme !----------------------------------------------------------------------- ktopcd = 0 kendcd = 0 tracer_dt = 0. wdep_in = 0. wdep_bc = 0. beta = 0. reevap_fraction = 0. reevap_diag = 0. tracer_flux = 0. sum_wdep = 0. id = size(T,1) jd = size(T,2) kd = size(T,3) call get_tracer_names(MODEL_ATMOS,n,tracer_name, units = units) if ( .not. Wetdep(n)%Lwetdep) return text_in_scheme = Wetdep(n)%text_in_scheme control = Wetdep(n)%control scheme = Wetdep(n)%scheme Henry_constant = Wetdep(n)%Henry_constant Henry_variable = Wetdep(n)%Henry_variable frac_in_cloud = Wetdep(n)%frac_in_cloud frac_in_cloud_snow = Wetdep(n)%frac_in_cloud_snow alpha_r = Wetdep(n)%alpha_r alpha_s = Wetdep(n)%alpha_s Lwetdep = Wetdep(n)%Lwetdep Lgas = Wetdep(n)%Lgas Laerosol = Wetdep(n)%Laerosol Lice = Wetdep(n)%Lice rho_air(:,:,:) = pfull(:,:,:) / ( T(:,:,:)*RDGAS ) ! kg/m3 ! Lice = .not. (scheme=='henry_noice' .or. scheme=='henry_below_noice' .or. & ! scheme=='aerosol_noice' .or. scheme=='aerosol_below_noice' ) if (Lice) then rainsnow3d(:,:,:) = rain3d(:,:,:) + snow3d(:,:,:) else rainsnow3d(:,:,:) = rain3d(:,:,:) end if do k=1,kd precip3d(:,:,k) = rainsnow3d(:,:,k+1)-rainsnow3d(:,:,k) precip3dr(:,:,k) = rain3d(:,:,k+1)-rain3d(:,:,k) pwt(:,:,k) = ( phalf(:,:,k+1) - phalf(:,:,k) )/GRAV ! kg/m2 zdel(:,:,k) = zhalf(:,:,k) - zhalf(:,:,k+1) ! m end do if (Lice) then do k=1,kd precip3ds(:,:,k) = snow3d(:,:,k+1)-snow3d(:,:,k) end do else do k=1,kd precip3ds(:,:,k) = 0. end do endif ! !+++ pag: ! if(lowercase(scheme)=='wdep_gas' .or. lowercase(scheme)=='wdep_aerosol') then ph = 5.0 ! cloud liquid water content xliq(:,:,:) = 0.5E-3 !default value if(trim(cloud_param) .eq. 'convect') then xliq(:,:,:) = 1.0e-3 elseif(trim(cloud_param) .eq. 'lscale') then xliq(:,:,:) = 0.5e-3 endif ! ! tracer == gas ! if(lowercase(scheme)=="wdep_gas") THEN frac_int = frac_int_gas do k = 1, kd do j = 1, jd do i = 1, id Htemp(i,j,k)=henry_constant & *exp(-Henry_variable*(1./298.-1./T(i,j,k))) K1(i,j)=1.2e-2*exp(-2010*(1/298.-1/T(i,j,k))) K2(i,j)=6.6e-8*exp(-1510*(1/298.-1/T(i,j,k))) HTemp(i,j,k)=Htemp(i,j,k)*(1 + K1(i,j) & /10.**(-ph) + K1(i,j)*K2(i,j)/(10.**(-ph))**2) f_a(i,j) = Htemp(i,j,k)/101.325*RGAS & *T(i,j,k)*xliq(i,j,k)*rho_air(i,j,k)/DENS_H2O scav_fact3d(i,j,k)=f_a(i,j)/(1.+f_a(i,j)) enddo enddo enddo elseif(lowercase(scheme)=="wdep_aerosol") THEN frac_int = frac_int_aerosol scav_fact3d(:,:,:)=frac_in_cloud else print *,' Aerosol number =',n,' tracer_name=',tracer_name,' scheme=',text_in_scheme print *, 'Please check "am2p12.ft' call ERROR_MESG('wet_deposition', 'Tracer is neither aerosol NOR gas.', FATAL ) endif ! !in cloud scavenging ! do k=1,kd do j = 1, jd do i = 1, id beta(i,j) = MAX( 0.0, precip3d(i,j,k)/pwt(i,j,k)/xliq(i,j,k)) in_temp(i,j,k) = (exp(-beta(i,j)*scav_fact3d(i,j,k)*dt)-1.0) if ( tracer(i,j,k) .gt. 0.) then wdep_in(i,j)=wdep_in(i,j) & - in_temp(i,j,k)*tracer(i,j,k)*pwt(i,j,k) tracer_dt(i,j,k) = tracer_dt(i,j,k) & - in_temp(i,j,k)*tracer(i,j,k)/dt endif ! !--reevaporation ! calculation of fracion of aerosols to-be ! reevaporated to the atmosphere: beta(i,j)=precip3d(i,j,k) if (beta(i,j) < 0.) then beta(i,j) = beta(i,j)/rainsnow3d(i,j,k) endif if (rainsnow3d(i,j,k+1) == 0. ) then !--reevaporation total beta(i,j)=MIN(MAX(0.,-beta(i,j)),1.) else beta(i,j)=MIN(MAX(0.,-beta(i,j))*frac_int,1.) endif ! reevporating to atmosphere reevap_diag(i,j,k)=beta(i,j)*wdep_in(i,j) wdep_in(i,j) = wdep_in(i,j)*(1.-beta(i,j)) tracer_dt(i,j,k) = tracer_dt(i,j,k) & - reevap_diag(i,j,k)/pwt(i,j,k)/dt enddo enddo enddo ! Below cloud scavenging do k=1,kd do j = 1, jd do i = 1, id fluxs(i,j) = (snow3d(i,j,k+1)+snow3d(i,j,k))/2.0 fluxr(i,j) = (rain3d(i,j,k+1)+rain3d(i,j,k))/2.0 bc_temp(i,j,k) = 3./4.*dt* & (fluxr(i,j)*alpha_r/R_r/DENS_H2O & + fluxs(i,j)*alpha_s/R_s/DENS_SNOW) if ( tracer(i,j,k) .gt. 0. ) then wdep_bc(i,j)=wdep_bc(i,j) & + bc_temp(i,j,k)*tracer(i,j,k)*pwt(i,j,k) tracer_dt(i,j,k) = tracer_dt(i,j,k) & + bc_temp(i,j,k)*tracer(i,j,k)/dt endif enddo enddo enddo ! ! end wdep_gas or wdep_aerosol ! else ! Calculate fraction of precipitation reevaporated in layer do k=1,kd where( rainsnow3d(:,:,k) > 0. .and. precip3d(:,:,k) < 0. ) reevap_fraction(:,:,k) = & -precip3d(:,:,k) / (rainsnow3d(:,:,k)) ! fraction end where ! Assume that the tracer reevaporation fraction is 50% of the precip ! reevaporation fraction, except when fraction = 100% ! where( reevap_fraction(:,:,k) < 1. ) ! reevap_fraction(:,:,k) = 0.5*reevap_fraction(:,:,k) ! end where end do ! cloud liquid water content ! xliq = 0.5E-3 !default value ! if(trim(cloud_param) .eq. 'convect') then ! xliq = 1.0e-3 ! elseif(trim(cloud_param) .eq. 'lscale') then ! xliq = 0.5e-3 ! endif ! Lgas = lowercase(scheme)=='henry' .or. lowercase(scheme)=='henry_below' .or. & ! lowercase(scheme)=='henry_noice' .or. lowercase(scheme)=='henry_below_noice' ! Laerosol = lowercase(scheme)=='aerosol' .or. lowercase(scheme)=='aerosol_below' .or. & ! lowercase(scheme)=='aerosol_noice' .or. lowercase(scheme)=='aerosol_below_noice' ! Assume that the aerosol reevaporation fraction is 50% of the precip ! reevaporation fraction, except when fraction = 100% if( Lgas ) then frac_int = frac_int_gas elseif( Laerosol ) then frac_int = frac_int_aerosol else frac_int = 1. end if if( Lgas .or. Laerosol ) then ! units = VMR ! ! Henry_constant (mole/L/Pa) = [X](aq) / Px(g) ! where [X](aq) is the concentration of tracer X in precipitation (mole/L) ! Px(g) is the partial pressure of the tracer in the air (Pa) ! ! VMR (total) = VMR (gas) + VMR (aq) ! = VMR (gas) + [X] * L ! ! where L = cloud liquid amount (kg H2O/mole air) ! ! Using Henry's Law, [X] = H * Px = H * VMR(gas) * Pfull ! ! So, VMR (total) = VMR(gas) * [ 1 + H * Pfull * L ] ! ! VMR(gas) = VMR(total) / [1 + H * Pfull * L] ! ! [X] = H * Pfull * VMR(total) / [ 1 + H * Pfull * L] ! ! Following Giorgi and Chameides, JGR, 90(D5), 1985, the first-order loss ! rate constant (s^-1) of X due to wet deposition equals: ! ! k = W_X / n_X ! ! where W_x = the loss rate (molec/cm3/s), and n_X = the number density (molec/cm3) ! ! W_X = [X] * W_H2O / (55 mole/L) ! n_x = VMR(total) * n_air (molec/cm3) = VMR(total) * P/(kT) * 1E-6 m3/cm3 ! ! where P = atmospheric pressure (Pa) ! k = Boltzmann's constant = 1.38E-23 J/K ! T = temperature (K) ! W_H2O = removal rate of water (molec/cm3/s) ! ! [X] * W_H2O / 55 ! So, k = ------------------------------ ! VMR(total) * P/(kT) * 1E-6 ! ! W_H2O H * VMR(total) * P / [ 1 + H * P *L ] ! = ----- * --------------------------------------- ! 55 VMR(total) * P/(kT) * 1E-6 ! ! W_H2O H * kT * 1E6 ! = ----- * ------------- ! 55 1 + H * P * L ! ! W_H2O 1 1 H * P * L ! = ----- * ----- * - * ------------- ! 55 n_air L 1 + H * P * L ! ! where W_H2O = precip3d (kg/m2/s) * (AVOGNO/mw_h2o) (molec/kg) / zdel (m) * 1E-6 m3/cm3 ! if( (Lgas .and. Henry_constant > 0) .or. Laerosol ) then in_temp(:,:,:) = 0. bc_temp(:,:,:) = 0. do k=1,kd ! Calculate the temperature dependent Henry's Law constant scav_factor(:,:) = 0.0 xliq(:,:,k) = MAX( cloud(:,:,k) * mw_air, 0. ) ! (kg H2O)/(mole air) n_air(:,:,k) = pfull(:,:,k) / (kboltz*T(:,:,k)) * cm3_2_m3 ! molec/cm3 if (Lgas) then temp_factor(:,:) = 1/T(:,:,k)-inv298p15 Htemp(:,:,k) = Henry_constant * & exp( Henry_variable*temp_factor ) f_a(:,:) = Htemp(:,:,k) * pfull(:,:,k) * xliq(:,:,k) ! / cloud_frac scav_factor(:,:) = f_a(:,:) / ( 1.+f_a(:,:) ) scav_factor_s(:,:) = scav_factor(:,:) else if (Laerosol) then scav_factor(:,:) = frac_in_cloud if (frac_in_cloud == frac_in_cloud_snow) then scav_factor_s(:,:) = scav_factor(:,:) else scav_factor_s(:,:) = f_snow_berg(:,:,k)*frac_in_cloud_snow +& (1.-f_snow_berg(:,:,k))*frac_in_cloud endif end if ! where (precip3d(:,:,k) > 0.0) where (precip3d(:,:,k) > 0. .and. xliq(:,:,k) > 0.) w_h2o(:,:) = precip3d(:,:,k) * (AVOGNO/mw_h2o) / zdel(:,:,k) * cm3_2_m3 ! molec/cm3/s beta(:,:) = w_h2o(:,:) * mw_h2o / (n_air(:,:,k) * xliq(:,:,k)) where (precip3ds(:,:,k) > 0.0 .and. precip3dr(:,:,k) > 0.0) where (scav_factor(:,:) /= scav_factor_s(:,:)) scav_factor(:,:) = ( scav_factor(:,:)*precip3dr(:,:,k) + & scav_factor_s(:,:) * precip3ds(:,:,k)) / precip3d(:,:,k) end where elsewhere where (precip3ds(:,:,k) > 0.0) scav_factor(:,:) = scav_factor_s(:,:) endwhere endwhere in_temp(:,:,k) = beta(:,:) * scav_factor(:,:) ! 1/s endwhere enddo !----------------------------------------------------------------- ! Below-cloud wet scavenging !----------------------------------------------------------------- if( lowercase(scheme)=='henry_below' .or. lowercase(scheme)=='henry_below_noice') then k_g = d_g/rain_diam * & ( 2. + 0.6 * sqrt( rain_diam*rain_vterm/vk_air ) * (vk_air/d_g)**(1./3.) ) do i = 1,id do j = 1,jd conc(:) = tracer(i,j,:) * n_air(i,j,:) / cm3_2_m3 ! Convert from VMR to molec/m3 do kk = 1,kd stay = 1. if( precip3d(i,j,kk) > 0. ) then conc_rain_total = 0. stay = zfull(i,j,kk) / (rain_vterm * dt) stay = min( stay, 1. ) do k = kk,kd f_a0 = Htemp(i,j,k) * pfull(i,j,k) * xliq(i,j,kk) * n_air(i,j,kk)/n_air(i,j,k) scav_factor0 = f_a0 / ( 1.+f_a0 ) conc_sat = conc(k) * scav_factor0 ! molec/m3 <== (xeqca1) sa_drop0 = geo_fac / rain_diam * xliq(i,j,kk) * n_air(i,j,kk) / & ( DENS_H2O * AVOGNO * cm3_2_m3 ) ! (m2 H2O) / (m3 air) fgas0 = conc(k) * k_g ! molec/m2/s fall_time = zdel(i,j,k) / rain_vterm ! sec conc_rain = fgas0 * sa_drop0 * fall_time ! molec/m3 <== (xca1) conc_rain_total = conc_rain_total + conc_rain ! molec/m3 <== (all1) if ( conc_rain_total < conc_sat ) then conc(k) = max( conc(k)-conc_rain, 0. ) end if end do conc(kk) = conc(kk) / n_air(i,j,kk) * cm3_2_m3 ! Convert to VMR conc(kk) = tracer(i,j,kk) - conc(kk) if ( conc(kk) /= 0. .and. tracer(i,j,kk) /= 0. ) then fall_time = zdel(i,j,kk)/rain_vterm bc_temp(i,j,kk) = bc_temp(i,j,kk) + & conc(kk) / (tracer(i,j,kk) * fall_time) * stay ! 1/s end if end if end do end do end do else if ( lowercase(scheme) == 'aerosol_below' .or. lowercase(scheme) == 'aerosol_below_noice') then do k=1,kd fluxs = (snow3d(:,:,k+1)+snow3d(:,:,k))/2.0 fluxr = (rain3d(:,:,k+1)+rain3d(:,:,k))/2.0 bc_temp(:,:,k) = 3./4. * & (fluxr(:,:)*alpha_r/R_r/DENS_H2O + & fluxs(:,:)*alpha_s/R_s/DENS_SNOW) end do end if do k = 1,kd wdep_in(:,:) = wdep_in(:,:) - & in_temp(:,:,k)*tracer(:,:,k)*pwt(:,:,k)* & cloud_frac(:,:,k) wdep_bc(:,:) = wdep_bc(:,:) - & bc_temp(:,:,k)*tracer(:,:,k)*pwt(:,:,k) enddo dt_temp(:,:,:) = 1. - exp( -bc_temp(:,:,:)*dt ) & ! fractional loss/timestep * ( cloud_frac(:,:,:)*exp( -in_temp(:,:,:)*dt ) + (1-cloud_frac(:,:,:)) ) tracer_dt(:,:,:) = dt_temp(:,:,:) / dt !+ve loss frequency (1/sec) endif else if(lowercase(scheme)=='fraction') then tracer_dt = 0.0 !----------------------------------------------------------------------- ! ! Compute areal fractions experiencing wet deposition: ! ! Set minimum precipitation rate below which no wet removal ! occurs to 0.01 cm/day ie 1.16e-6 mm/sec (kg/m2/s) premin=1.16e-6 ! ! Large scale cloud liquid water content (kg/m3) ! and below cloud washout efficiency (cm-1): flaglw =parse(control,'lslwc',clwc) if (flaglw == 0 ) clwc=0.5e-3 wash=1.0 ! ! When convective adjustment occurs, use convective cloud liquid water content: ! if(trim(cloud_param) .eq. 'convect') then flaglw = parse(control,'convlwc',clwc) if (flaglw == 0) clwc=2.0e-3 wash=0.3 end if ! do j=1,size(rain,2) do i=1,size(rain,1) tracer_dt(i,j,:)=0.0 washout(i,j)=0.0 prenow = rain(i,j) + snow(i,j) if(prenow .gt. premin) then ! ! Assume that the top of the cloud is where the highest model level ! specific humidity is reduced. And the the bottom of the cloud is the ! lowest model level where specific humidity is reduced. ! ktopcd(i,j) = 0 do k = kd,1,-1 if (qdt(i,j,k) < 0.0 ) ktopcd(i,j) = k enddo kendcd(i,j) = 0 do k = 1,kd if (qdt(i,j,k) < 0.0 ) kendcd(i,j) = k enddo ! ! Thickness of precipitating cloud deck: ! if(ktopcd(i,j).gt.1) then hwtop = 0.0 do k=ktopcd(i,j),kendcd(i,j) hwtop=hwtop+(phalf(i,j,k+1)-phalf(i,j,k))*rdgas*T(i,j,k)/grav/pfull(i,j,k) enddo do k=ktopcd(i,j),kendcd(i,j) ! Areal fraction affected by precip clouds (max = 0.5): tracer_dt(i,j,k)=prenow/(clwc*hwtop) end do endif washout(i,j)=prenow*wash endif end do end do endif ! Now multiply by the tracer mixing ratio to get the actual tendency. tracer_dt(:,:,:) = MIN( MAX(tracer_dt(:,:,:), 0.0E+00), 0.5/dt) where (tracer > 0.) tracer_dt = tracer_dt*tracer else where tracer_dt = 0. end where !++lwh ! ! Re-evaporation ! do k = 1,kd where (reevap_fraction(:,:,k) > 0.) reevap_diag(:,:,k) = reevap_fraction(:,:,k) * tracer_flux(:,:) ! tracer reevaporation fraction is reduced from precip reevaporation, ! except when complete reevaporation occurs where( reevap_fraction(:,:,k) < 1. ) reevap_diag(:,:,k) = reevap_diag(:,:,k) * frac_int end where tracer_dt(:,:,k) = tracer_dt(:,:,k) - reevap_diag(:,:,k) / pwt(:,:,k) end where tracer_flux(:,:) = tracer_flux(:,:) + tracer_dt(:,:,k)*pwt(:,:,k) end do !--lwh ! endif ! End branching pag/lwh ! ! Output diagnostics in kg/m2/s (if MMR) or mole/m2/s (if VMR) if(trim(units) .eq. 'mmr') then diag_scale = 1. elseif(trim(units) .eq. 'vmr') then diag_scale = mw_air ! kg/mole else write(*,*) ' Tracer number =',n,' tracer_name=',tracer_name write(*,*) ' scheme=',text_in_scheme write(*,*) ' control=',control write(*,*) ' scheme=',scheme write(*,*) 'Please check field table' write(*,*) 'tracers units =',trim(units),'it should be either mmr or vmr!' ! ! Tracer units must be either VMR or MMR ! call error_mesg('wet_deposition', 'Unsupported tracer units.', FATAL ) endif ! Column integral of wet deposition sum_wdep = 0. do k=1,kd sum_wdep = sum_wdep + tracer_dt(:,:,k)*pwt(:,:,k)/diag_scale end do if (present (sum_wdep_out)) sum_wdep_out = -sum_wdep if(trim(cloud_param) == 'lscale') then if (id_tracer_reevap_ls(n) > 0 ) then used = send_data ( id_tracer_reevap_ls(n), reevap_diag/diag_scale, Time ,is,js,1) endif if (id_tracer_wdep_ls(n) > 0 ) then used = send_data ( id_tracer_wdep_ls(n), sum_wdep, Time, is_in =is, js_in=js ) endif if (id_tracer_wdep_ls_3d(n) > 0 ) then used = send_data ( id_tracer_wdep_ls_3d(n), tracer_dt*pwt/diag_scale, Time, is_in =is, js_in=js ,ks_in=1) endif if (id_tracer_wdep_lsin(n) > 0 ) then used = send_data ( id_tracer_wdep_lsin(n), wdep_in/diag_scale, Time, is_in=is, js_in=js ) endif if (id_tracer_wdep_lsbc(n) > 0 ) then used = send_data ( id_tracer_wdep_lsbc(n), wdep_bc/diag_scale, Time, is_in=is, js_in=js ) endif else if(trim(cloud_param) == 'convect') then if (id_tracer_reevap_cv(n) > 0 ) then used = send_data ( id_tracer_reevap_cv(n), reevap_diag/diag_scale, Time ,is,js,1) endif if(id_tracer_wdep_cv(n) > 0) then used = send_data( id_tracer_wdep_cv(n), sum_wdep, Time, is_in=is, js_in=js) endif if (id_tracer_wdep_cvin(n) > 0 ) then used = send_data ( id_tracer_wdep_cvin(n), wdep_in/diag_scale, Time, is_in=is, js_in=js) endif if (id_tracer_wdep_cvbc(n) > 0 ) then used = send_data ( id_tracer_wdep_cvbc(n), wdep_bc/diag_scale, Time, is_in=is, js_in=js) endif endif end subroutine wet_deposition ! ! !####################################################################### ! subroutine get_drydep_param(text_in_scheme,text_in_param,scheme,land_dry_dep_vel,sea_dry_dep_vel) !subroutine get_drydep_param(text_in_scheme, text_in_param, scheme, & ! land_dry_dep_vel, sea_dry_dep_vel, & ! ice_dry_dep_vel, snow_dry_dep_vel, & ! vegn_dry_dep_vel) ! ! Subroutine to initialiize the parameters for the dry deposition scheme. ! If the dry dep scheme is a "fixed" form then the ! dry_deposition velocity value has to be set. ! If the dry dep scheme is a "wind_driven" form then the dry_deposition ! velocity value will be calculated. So set to a dummy value of 0.0 ! INTENT IN ! text_in_scheme : The text that has been parsed from tracer table as ! the dry deposition scheme to be used. ! text_in_param : The parameters that are associated with the dry ! deposition scheme. ! INTENT OUT ! scheme : The scheme that is being used. ! land_dry_dep_vel : Dry deposition velocity over the land ! sea_dry_dep_vel : Dry deposition velocity over the sea ! character(len=*), intent(in) :: text_in_scheme, text_in_param character(len=*), intent(out) :: scheme real, intent(out) :: land_dry_dep_vel, sea_dry_dep_vel!, & ! ice_dry_dep_vel, snow_dry_dep_vel, & ! vegn_dry_dep_vel integer :: flag !Default scheme = 'None' land_dry_dep_vel=0.0 sea_dry_dep_vel=0.0 !ice_dry_dep_vel=0.0 !snow_dry_dep_vel=0.0 !vegn_dry_dep_vel=0.0 if(lowercase(trim(text_in_scheme(1:4))).eq.'wind') then scheme = 'Wind_driven' endif !if(lowercase(trim(text_in_scheme(1:15))).eq.'sfc_dependent_w') then !scheme = 'sfc_dependent_wind_driven' !endif ! !if(lowercase(trim(text_in_scheme(1:6))).eq.'sfc_bl') then !scheme = 'sfc_BL_dependent_wind_driven' !endif if(lowercase(trim(text_in_scheme(1:5))).eq.'fixed') then scheme = 'fixed' flag=parse(text_in_param,'land',land_dry_dep_vel) flag=parse(text_in_param,'sea', sea_dry_dep_vel) endif !if(lowercase(trim(text_in_scheme(1:15))).eq.'sfc_dependent_f') then !scheme = 'sfc_dependent_fixed' !flag=parse(text_in_param,'land',land_dry_dep_vel) !flag=parse(text_in_param,'sea', sea_dry_dep_vel) !flag=parse(text_in_param,'ice', ice_dry_dep_vel) !if (flag == 0) then ! ice_dry_dep_vel = sea_dry_dep_vel !end if !flag=parse(text_in_param,'snow', snow_dry_dep_vel) !if (flag == 0) then ! snow_dry_dep_vel = land_dry_dep_vel !end if !flag=parse(text_in_param,'vegn', vegn_dry_dep_vel) !if (flag == 0) then ! vegn_dry_dep_vel = land_dry_dep_vel !end if !endif if(lowercase(trim(text_in_scheme(1:4))).eq.'file') then scheme = 'file' endif end subroutine get_drydep_param ! !####################################################################### ! ! ! subroutine get_wetdep_param(text_in_scheme,text_in_param,scheme,& henry_constant, henry_temp, & frac_in_cloud, frac_in_cloud_snow, & alpha_r,alpha_s, & Lwetdep, Lgas, Laerosol, Lice, & frac_in_cloud_uw, frac_in_cloud_donner) ! ! Routine to initialize the parameters for the wet deposition scheme. ! ! ! shm has modified this subroutine to include additional parameters: ! frac_in_cloud, alpha_r, and alpha_s ! ! INTENT IN ! ! Text read from the tracer table which provides information on which ! wet deposition scheme to use. ! ! ! Parameters associated with the wet deposition scheme. These will be ! parsed in this routine. ! ! ! Wet deposition scheme to use. ! Choices are: None, Fraction, Henry, Henry_below, Aerosol, Aerosol_below, ! wdep_aerosol, wdep_gas ! ! ! Henry's Law constant for the tracer (see wet_deposition for explanation of Henry's Law) ! ! ! The temperature dependence of the Henry's Law constant. ! ! ! In-cloud fraction for aerosols ! ! ! Controls below-cloud aerosol scavenging by rain ! ! ! Controls below-cloud aerosol scavenging by snow ! ! ! Does tracer have wet removal? ! ! ! Is tracer a gas? ! ! ! Is tracer an aerosol? ! ! ! Is tracer removed by snow (or just rain)? ! character(len=*), intent(in) :: text_in_scheme, text_in_param character(len=*), intent(out) :: scheme real, intent(out) :: henry_constant, henry_temp real, intent(out) :: frac_in_cloud, frac_in_cloud_snow real, intent(out) :: alpha_r, alpha_s logical, intent(out) :: Lwetdep, Lgas, Laerosol, Lice real, intent(out), optional :: frac_in_cloud_uw, frac_in_cloud_donner integer :: flag !Default scheme = 'None' henry_constant= 0. henry_temp = 0. frac_in_cloud = 0. frac_in_cloud_snow = 0. alpha_r = 0. alpha_s = 0. Lwetdep = .false. Lgas = .false. Laerosol = .false. if (present(frac_in_cloud_uw)) frac_in_cloud_uw = 0. if (present(frac_in_cloud_donner)) frac_in_cloud_donner = 0. if( trim(lowercase(text_in_scheme)) == 'fraction' ) then scheme = 'Fraction' else if( trim(lowercase(text_in_scheme)) == 'henry' .or. & trim(lowercase(text_in_scheme)) == 'henry_below' .or. & trim(lowercase(text_in_scheme)) == 'henry_noice' .or. & trim(lowercase(text_in_scheme)) == 'henry_below_noice' ) then if( trim(lowercase(text_in_scheme)) == 'henry' ) then scheme = 'henry' else if ( trim(lowercase(text_in_scheme)) == 'henry_below' ) then scheme = 'henry_below' else if ( trim(lowercase(text_in_scheme)) == 'henry_noice' ) then scheme = 'henry_noice' else if ( trim(lowercase(text_in_scheme)) == 'henry_below_noice' ) then scheme = 'henry_below_noice' end if flag=parse(text_in_param,'henry', henry_constant) flag=parse(text_in_param,'dependence',henry_temp ) Lgas = .true. else if( trim(lowercase(text_in_scheme)) == 'aerosol' .or. & trim(lowercase(text_in_scheme)) == 'aerosol_below' .or. & trim(lowercase(text_in_scheme)) == 'aerosol_noice' .or. & trim(lowercase(text_in_scheme)) == 'aerosol_below_noice' ) then if( trim(lowercase(text_in_scheme)) == 'aerosol' ) then scheme = 'aerosol' else if ( trim(lowercase(text_in_scheme)) == 'aerosol_below' ) then scheme = 'aerosol_below' else if ( trim(lowercase(text_in_scheme)) == 'aerosol_noice' ) then scheme = 'aerosol_noice' else if ( trim(lowercase(text_in_scheme)) == 'aerosol_below_noice' ) then scheme = 'aerosol_below_noice' end if flag=parse(text_in_param,'frac_incloud',frac_in_cloud) flag=parse(text_in_param,'frac_incloud_snow',frac_in_cloud_snow) if (flag == 0) then frac_in_cloud_snow = frac_in_cloud end if if (present(frac_in_cloud_uw)) then flag=parse(text_in_param,'frac_incloud_uw',frac_in_cloud_uw) if (flag == 0) then frac_in_cloud_uw = frac_in_cloud end if end if if (present(frac_in_cloud_donner)) then flag=parse(text_in_param,'frac_incloud_donner', & frac_in_cloud_donner) if (flag == 0) then frac_in_cloud_donner = frac_in_cloud end if end if flag=parse(text_in_param,'alphar',alpha_r) flag=parse(text_in_param,'alphas',alpha_s) Laerosol = .true. end if if( trim(lowercase(text_in_scheme)) == 'wdep_aerosol') scheme= 'wdep_aerosol' if ( trim(lowercase(text_in_scheme)) == 'wdep_gas' ) scheme= 'wdep_gas' if (scheme .eq. 'wdep_aerosol' .or. scheme .eq. 'wdep_gas') then flag=parse(text_in_param,'frac_incloud',frac_in_cloud) flag=parse(text_in_param,'alphar',alpha_r) flag=parse(text_in_param,'alphas',alpha_s) end if Lice = .not. ( scheme=='henry_noice' .or. scheme=='henry_below_noice' .or. & scheme=='aerosol_noice' .or. scheme=='aerosol_below_noice' ) Lwetdep = scheme /= 'None' end subroutine get_wetdep_param ! ! !####################################################################### ! ! subroutine interp_emiss(global_source, start_lon, start_lat, & lon_resol, lat_resol, data_out) ! ! ! A routine to interpolate emission fields of arbitrary resolution onto the ! resolution of the model. ! ! ! Routine to interpolate emission fields (or any 2D field) to the model ! resolution. The local section of the global field is returned to the ! local processor. ! ! ! ! INTENT IN ! ! Global emission field. ! ! ! Longitude of starting point of emission field ! (in radians). This is the westernmost boundary of the ! global field. ! ! ! Latitude of starting point of emission field ! (in radians). This is the southern boundary of the ! global field. ! ! ! Longitudinal resolution of the emission data (in radians). ! ! ! Latitudinal resolution of the emission data (in radians). ! ! ! INTENT OUT ! ! Interpolated emission field on the local PE. ! real, intent(in) :: global_source(:,:) real, intent(in) :: start_lon,start_lat,lon_resol,lat_resol real, intent(out) :: data_out(:,:) integer :: i, j, nlon_in, nlat_in real :: blon_in(size(global_source,1)+1) real :: blat_in(size(global_source,2)+1) type (horiz_interp_type) :: Interp ! Set up the global surface boundary condition longitude-latitude boundary values nlon_in = size(global_source,1) nlat_in = size(global_source,2) ! For some reason the input longitude needs to be incremented by 180 degrees. do i = 1, nlon_in+1 blon_in(i) = start_lon + float(i-1)*lon_resol + PI enddo if (abs(blon_in(nlon_in+1)-blon_in(1)-twopi) < epsilon(blon_in)) & blon_in(nlon_in+1)=blon_in(1)+twopi do j = 2, nlat_in blat_in(j) = start_lat + float(j-1)*lat_resol enddo blat_in(1) = -0.5*PI blat_in(nlat_in+1) = 0.5*PI ! Now interpolate the global data to the model resolution call horiz_interp_init call horiz_interp_new (Interp, blon_in, blat_in, & blon_out, blat_out) call horiz_interp (Interp, global_source, data_out) call horiz_interp_del ( Interp ) end subroutine interp_emiss ! ! ###################################################################### ! subroutine GET_RH (T,Q,P,RH,mask) !*********************************************************************** ! SUBROUTINE GET_RH ! PURPOSE ! VECTOR COMPUTATION OF RELATIVE HUMIDITY ! DESCRIPTION OF PARAMETERS ! T TEMPERATURE VECTOR (DEG K) ! P PRESSURE VECTOR (Pa) ! Q SPECIFIC HUMIDITY (kg/kg) ! RH RELATIVE HUMIDITY ! MASK EARTH SURFACE BELOW GROUND (in pressure coordinates) ! !*********************************************************************** ! Real, parameter :: ONE = 1. Real, parameter :: ZP622 = 0.622 Real, parameter :: Z1P0S1 = 1.00001 Real, parameter :: Z1P622 = 1.622 Real, parameter :: Z138P9 = 138.90001 Real, parameter :: Z198P9 = 198.99999 Real, parameter :: Z200 = 200.0 Real, parameter :: Z337P9 = 337.9 real, intent(in), dimension(:,:,:) :: T !temp at curr time step [ deg k ] real, intent(in), dimension(:,:,:) :: P !pressure at full levels [ Pa ] real, intent(in), dimension(:,:,:) :: Q !specific humidity at current time step [ kg / kg ] real, intent(in), dimension(:,:,:), optional :: mask ! real, intent(out), dimension(:,:,:) :: RH !relative humidity [0-1] ! Dynamic Work Space ! ------------------ real :: A1622 real, dimension(size(q,1),size(q,2),size(q,3)) :: e1, e2, tq, qs integer, dimension(size(q,1),size(q,2),size(q,3)) :: i1, i2 integer :: i,j,k ! real, dimension(67) :: EST1 data EST1/ 0.31195E-02, 0.36135E-02, 0.41800E-02, & 0.48227E-02, 0.55571E-02, 0.63934E-02, 0.73433E-02, & 0.84286E-02, 0.96407E-02, 0.11014E-01, 0.12582E-01, & 0.14353E-01, 0.16341E-01, 0.18574E-01, 0.21095E-01, & 0.23926E-01, 0.27096E-01, 0.30652E-01, 0.34629E-01, & 0.39073E-01, 0.44028E-01, 0.49546E-01, 0.55691E-01, & 0.62508E-01, 0.70077E-01, 0.78700E-01, 0.88128E-01, & 0.98477E-01, 0.10983E+00, 0.12233E+00, 0.13608E+00, & 0.15121E+00, 0.16784E+00, 0.18615E+00, 0.20627E+00, & 0.22837E+00, 0.25263E+00, 0.27923E+00, 0.30838E+00, & 0.34030E+00, 0.37520E+00, 0.41334E+00, 0.45497E+00, & 0.50037E+00, 0.54984E+00, 0.60369E+00, 0.66225E+00, & 0.72589E+00, 0.79497E+00, 0.86991E+00, 0.95113E+00, & 0.10391E+01, 0.11343E+01, 0.12372E+01, 0.13484E+01, & 0.14684E+01, 0.15979E+01, 0.17375E+01, 0.18879E+01, & 0.20499E+01, 0.22241E+01, 0.24113E+01, 0.26126E+01, & 0.28286E+01, 0.30604E+01, 0.33091E+01, 0.35755E+01/ ! real, dimension(72) :: EST2 data EST2/ & 0.38608E+01, 0.41663E+01, 0.44930E+01, 0.48423E+01, & 0.52155E+01, 0.56140E+01, 0.60394E+01, 0.64930E+01, & 0.69767E+01, 0.74919E+01, 0.80406E+01, 0.86246E+01, & 0.92457E+01, 0.99061E+01, 0.10608E+02, 0.11353E+02, & 0.12144E+02, 0.12983E+02, 0.13873E+02, 0.14816E+02, & 0.15815E+02, 0.16872E+02, 0.17992E+02, 0.19176E+02, & 0.20428E+02, 0.21750E+02, 0.23148E+02, 0.24623E+02, & 0.26180E+02, 0.27822E+02, 0.29553E+02, 0.31378E+02, & 0.33300E+02, 0.35324E+02, 0.37454E+02, 0.39696E+02, & 0.42053E+02, 0.44531E+02, 0.47134E+02, 0.49869E+02, & 0.52741E+02, 0.55754E+02, 0.58916E+02, 0.62232E+02, & 0.65708E+02, 0.69351E+02, 0.73168E+02, 0.77164E+02, & 0.81348E+02, 0.85725E+02, 0.90305E+02, 0.95094E+02, & 0.10010E+03, 0.10533E+03, 0.11080E+03, 0.11650E+03, & 0.12246E+03, 0.12868E+03, 0.13517E+03, 0.14193E+03, & 0.14899E+03, 0.15634E+03, 0.16400E+03, 0.17199E+03, & 0.18030E+03, 0.18895E+03, 0.19796E+03, 0.20733E+03, & 0.21708E+03, 0.22722E+03, 0.23776E+03, 0.24871E+03/ ! real, dimension(139) :: EST EQUIVALENCE (EST(1) , EST1(1)), (EST(68),EST2(1)) !*********************************************************************** ! A1622 = ONE / Z1P622 TQ = T - Z198P9 I1(:,:,:) = 1 I2(:,:,:) = 1 where ( T < Z200 ) TQ = Z1P0S1 where ( T > Z337P9 ) TQ = Z138P9 IF (present(mask)) THEN where ( mask > 0. ) I1 = int(TQ) I2 = I1 + 1 end where else I1 = int(TQ) I2 = I1 + 1 endif do i=1,size(q,1) do j=1,size(q,2) do k=1,size(q,3) E1(i,j,k) = EST( I1(i,j,k) ) E2(i,j,k) = EST( I2(i,j,k) ) enddo enddo enddo QS(:,:,:) = TQ(:,:,:) - float(I1(:,:,:)) QS(:,:,:) = E1(:,:,:) + QS(:,:,:) * ( E2(:,:,:)-E1(:,:,:) ) E1(:,:,:) = (0.01 * P(:,:,:)) * A1622 where ( E1 < QS ) QS = E1 if (present(mask)) then where ( mask > 0. ) QS = ZP622 * QS / ( P * 0.01) else QS(:,:,:) = ZP622 * QS(:,:,:) / ( P(:,:,:) * 0.01) endif RH(:,:,:) = Q(:,:,:)/QS(:,:,:) end subroutine GET_RH ! ###################################################################### ! subroutine get_w10m(z_full, u, v, rough_mom,u_star, b_star, q_star, & w10m_ocean, w10m_land, Time_next, is,js) real, intent(in), dimension(:,:) :: z_full, u, v real, intent(in), dimension(:,:) :: rough_mom real, intent(in), dimension(:,:) :: u_star, b_star, q_star type(time_type), intent(in) :: Time_next integer, intent(in) :: is,js logical :: used real, intent(out), dimension(:,:) :: w10m_ocean, w10m_land real, dimension(size(u,1),size(u,2)) :: del_m real, dimension(size(u,1),size(u,2)) :: del_h real, dimension(size(u,1),size(u,2)) :: del_q ! Reference heights for momentum and heat [m] real, parameter :: zrefm = 10. real, parameter :: zrefh = 2. real, parameter :: scaling_factor=1. w10m_ocean(:,:) = 0.0 w10m_land (:,:) = 0.0 del_m(:,:) = 0.0 call mo_profile(zrefm, zrefh, z_full, & rough_mom, rough_mom, rough_mom, & u_star, b_star, q_star, & del_m, del_h, del_q ) !----------------------------------------------------------------- ! ... Wind speed at anemometer level (10 meters above ground) !----------------------------------------------------------------- w10m_ocean(:,:)=sqrt(u(:,:)**2 +v(:,:)**2 )*del_m(:,:) w10m_land (:,:)=sqrt(u(:,:)**2 +v(:,:)**2 )*del_m(:,:)*scaling_factor ! Send the scaling factor if (id_delm > 0 ) then used = send_data ( id_delm, del_m, Time_next, is_in=is,js_in=js ) endif ! Send the 10m wind speed data to the diag_manager for output. if (id_w10m > 0 ) then used = send_data ( id_w10m, w10m_land, Time_next, is_in=is,js_in=js ) endif end subroutine get_w10m ! ###################################################################### ! subroutine get_cldf(ps, pfull, rh, cldf) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! !!! This subroutine estimates the cloud fraction "cldf" for !!! each grid box using an empirical function of the relative !!! humidity in that grid box, following Sundqvist et al., Mon. Weather Rev., !!! v117, 164101657, 1989: !!! !!! cldf = 1 - sqrt[ 1 - (RH - RH0)/(1 - RH0) ] !!! !!! where RH is the relative humidity and RH0 is the threshold relative !!! humidity for condensation specified as a function of pressure based !!! on Xu and Krueger, Mon. Weather Rev., v119, 342-367, 1991. !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real, intent(in), dimension(:,:) :: ps real, intent(in), dimension(:,:,:) :: pfull real, intent(in), dimension(:,:,:) :: rh real, intent(out), dimension(:,:,:) :: cldf real, parameter :: zrt = 0.6 real, parameter :: zrs = 0.99 integer :: i,j,k, id, jd, kd real :: p, r, r0, b0 id=size(pfull,1); jd=size(pfull,2); kd=size(pfull,3) do k = 1, kd do j = 1, jd do i = 1, id P = pfull(i,j,k) R = RH(i,j,k) R0 = ZRT + (ZRS-ZRT) * exp(1.-(PS(i,j)/P)**2.5) B0 = (R-R0) / (1.-R0) if (R .lt.R0) B0 = 0. if (B0.gt.1.) B0 = 1. CLDF(i,j,k) = 1.-sqrt(1.-B0) end do end do end do end subroutine get_cldf !###################################################################### ! ! ! The destructor routine for the tracer utilities module. ! ! ! This subroutine writes the version name to logfile and exits. ! subroutine atmos_tracer_utilities_end deallocate(blon_out, blat_out) module_is_initialized = .FALSE. end subroutine atmos_tracer_utilities_end ! ! ###################################################################### ! !IROUTINE: sjl_fillz --- Fill from neighbors below and above ! ! !INTERFACE: subroutine sjl_fillz(im, km, nq, q, dp) implicit none ! !INPUT PARAMETERS: integer, intent(in):: im ! No. of longitudes integer, intent(in):: km ! No. of levels integer, intent(in):: nq ! Total number of tracers real, intent(in):: dp(im,km) ! pressure thickness ! !INPUT/OUTPUT PARAMETERS: real, intent(inout) :: q(im,km,nq) ! tracer mixing ratio ! !DESCRIPTION: ! Check for "bad" data and fill from east and west neighbors ! ! !BUGS: ! Currently this routine only performs the east-west fill algorithm. ! This is because the N-S fill is very hard to do in a reproducible ! fashion when the problem is decomposed by latitudes. ! ! !REVISION HISTORY: ! 00.04.01 Lin Creation ! !EOP !----------------------------------------------------------------------- !BOC ! ! !LOCAL VARIABLES: integer i, k, ic real qup, qly, dup do ic=1,nq ! Top layer do i=1,im if( q(i,1,ic) < 0.) then q(i,2,ic) = q(i,2,ic) + q(i,1,ic)*dp(i,1)/dp(i,2) q(i,1,ic) = 0. endif enddo ! Interior do k=2,km-1 do i=1,im if( q(i,k,ic) < 0. ) then ! Borrow from above qup = q(i,k-1,ic)*dp(i,k-1) qly = -q(i,k ,ic)*dp(i,k ) dup = min( 0.75*qly, qup ) !borrow no more than 75% from top q(i,k-1,ic) = q(i,k-1,ic) - dup/dp(i,k-1) ! Borrow from below: q(i,k,ic) is still negative at this stage q(i,k+1,ic) = q(i,k+1,ic) + (dup-qly)/dp(i,k+1) q(i,k ,ic) = 0. endif enddo enddo ! Bottom layer k = km do i=1,im if( q(i,k,ic) < 0.) then ! Borrow from above qup = q(i,k-1,ic)*dp(i,k-1) qly = -q(i,k ,ic)*dp(i,k ) dup = min( qly, qup ) q(i,k-1,ic) = q(i,k-1,ic) - dup/dp(i,k-1) q(i,k,ic) = 0. endif enddo enddo end subroutine sjl_fillz end module atmos_tracer_utilities_mod