Update to silam_v5_7@596839

Bugfix in boundary-condition timing
Updated dust source

   git-svn-id: https://svn.fmi.fi/svn/tie/SILAM/silam_v5_7@596839 250f1e8d-2010-0410-9efa-e22dadc992bd

source/advection_eulerian_v5.silam.mod.f90

@@ -2760,7 +2760,7 @@ subroutine adv_diffusion_vertical_v5(pDispFlds, &
               ! Free bottom and top for diffusion       
               ! If we do not kill momens here -- we are in trouble
               passDiff(:,0) = 0.
-              passDiff(:,nz_dispersion+1) = 0.
+              passDiff(1:,nz_dispersion+1) = 0. !!! Keep mass above: Boundary-condition is there if diffusion goes first
 
 
                 if (diffuse_cm_vert)  &

source/dispersion_supplementary.silam.mod.f90

@@ -413,44 +413,13 @@ subroutine run_dispersion(cloud, em_source, wdr, simrules, &
             call stop_count(chCounterNm = command_string)
 
 
-!          bc_lim1 = fu_closest_obstime(reftime, backwards, fu_obstime_interval(simRules%IniBoundaryRules))
-!          bc_lim2 = fu_closest_obstime(reftime, forwards, fu_obstime_interval(simRules%IniBoundaryRules))
-!          IF(now == simrules%startTime)THEN 
-!            !--- For the first time data must always be acquired
-!            bc_lim1_prev = bc_lim1 + simRules%timestep
-!            bc_lim2_prev = bc_lim2 + simRules%timestep
-!          endif
-!          IF (error) EXIT loop_over_time
-!          IF(.not.(fu_between_times(bc_lim1, bc_lim1_prev, bc_lim2_prev, .true.) &
-!           & .and. fu_between_times(bc_lim2, bc_lim1_prev, bc_lim2_prev, .true.)))THEN
-!
-!            command_string = 'Boundary_condition_acquisition'
-!            call msg(command_string)
-!            call start_count(chCounterNm = command_string)
-!            !
-!            ! ATTENTION. There used to be now, now+timestep*2. in both time boundaries
-!            ! However, I have not found any reason to keep it.
-!            !
-!            do iTmp = 1, fu_nr_boundary_inFiles(simRules%IniBoundaryRules)
-!              CALL fix_shopping_time_boundaries(fu_shplst(simRules%IniBoundaryRules, iTmp), &
-!                                              & now, (now + simRules%timestep))
-!!call msg('')
-!!call msg('BC Shopping list:')
-!!call report(fu_shplst(simRules%IniBoundaryRules, iTmp))
-!!call msg('')
-!            enddo
-!            CALL fill_boundary_market(BCMarketPtr, simRules%IniBoundaryRules, first_step)
-!            IF (error) EXIT loop_over_time
-!
-!            !command_string = 'Boundary_condition_acquisition'
-!            call stop_count(chCounterNm = command_string)
 
             command_string = 'Boundary_condition_processing'
             call msg(command_string)
             call start_count(chCounterNm = command_string)
 
             CALL fillBoundaryStruct(simRules%IniBoundaryRules, BCMarketPtr, &
-                                  & now, fu_boundaryStructures(cloud), meteo_ptr, first_step)
+                                  & reftime, fu_boundaryStructures(cloud), meteo_ptr, first_step)
             IF (error) EXIT loop_over_time
 
             !command_string = 'Boundary_condition_processing'

source/ini_boundary_conditions.silam.mod.f90

@@ -137,7 +137,7 @@ MODULE ini_boundary_conditions
     real, dimension(:), allocatable :: polarCapThickness       ! (iVertical)
     real :: polarCapArea
     type(silja_field_3d_pointer), dimension(:,:), pointer :: Fld3dPtr   ! dimension - past/future; 1:nSpTr
-    integer :: past_future_switch = 0      ! 1 or 0                              
+    integer :: past_future_switch = 0      !  0 -> [past, future], 1 ->  [future -> past]
     type(THorizInterpStruct), pointer :: meteo2BoundaryInterpHoriz => null()
     type(TVertInterpStruct), pointer :: meteo2BoundaryInterpVert => null()
     type(silam_vertical) :: vertical
@@ -2065,7 +2065,7 @@ subroutine fillBoundaryStruct(ibRules, miniMarketBC, now, boundStructArray, &
       ! Imported parameters
       type(Tini_boundary_rules),  intent(inout), target :: ibRules
       type(TboundaryStructPtr), dimension(:), pointer :: boundStructArray
-      type(silja_time ), intent(in) :: now
+      type(silja_time ), intent(in) :: now !! Mid-time-step
       type(mini_market_of_stacks), pointer :: miniMarketBC
       type(Tfield_buffer), pointer :: met_buf
       logical, intent(in) :: ifFirstTime
@@ -2846,16 +2846,13 @@ subroutine get_dirichlet_boundary(bnd, buf_data, species_mapping_trn, nSpecies,
         integer, dimension(:), intent(in) :: species_mapping_trn
         integer, intent(in) :: x1, x2, y1, y2, z1, z2, nSpecies
         logical, intent(in) :: ifPlain
-        type(silja_time), intent(in) :: now
+        type(silja_time), intent(in) :: now  !!! Should be mid-step here
 
         ! Local variables
         integer :: ix, iy, iz, iCell, izb, iSpecies, iSpeciesTrn !, i1b
         real, dimension(max_species) :: weight_past_boundary
         real, dimension(:), pointer :: val_past, val_future !, &
-!                                     & z_cell_size_past, z_cell_size_future
-!        real :: weight_past
 
-!        weight_past = dusp_buf%weight_past
 
         !
         ! Get the weight_past_boundary, which can be species-specific

source/pollution_cloud.silam.mod.f90

@@ -1480,7 +1480,8 @@ subroutine advect_pollution_cloud_v4(cloud,now,timestep, &
       ! The main Eulerian advection
       !
       if(fu_ifBoundary(IniBoundaryRules))then
-        call boundary_conditions_now(cloud%boundaryStructures, cloud%pBoundaryBuffer, now)
+        !! Should be valid for mid-timestep
+        call boundary_conditions_now(cloud%boundaryStructures, cloud%pBoundaryBuffer, now + timestep*0.5)
         if(error)return
       endif
 #ifdef DEBUG

source/source_terms_wind_blown_dust.silam.mod.f90

@@ -129,7 +129,7 @@ MODULE source_terms_wind_blown_dust
 !                     & fSrfRoughSandOnly = 3.5e-5, &   ! too high emission over China, not sure about Africa. 
 !                     & fSrfRoughSandOnly = 4e-5, &   ! seems to be too high emission everywhere
                      & fDragPartitioningSandOnly = 1.0/log(0.35*(0.1/fSrfRoughSandOnly)**0.8), &
-                     & fLAI_critical = 1.0, &
+                     & fLAI_critical = 1.5, &
 !                     & u_star_excess_saturation = 0.1   ! m/s
 !                     & u_star_excess_saturation = 0.09   ! m/s
 !                     & u_star_excess_saturation = 0.06   ! m/s
@@ -953,7 +953,8 @@ subroutine project_wbdust_src_second_grd(wbdust_src, grid, vert_disp, vert_proj,
         arTmp(4) = 0.125
 
       case(simple_dust_flag)
-        call set_vertical(fu_set_layer_between_two(layer_btw_2_height, 1.0, 100.0), vertTmp)
+        ! All emission between 1 m and 500 m
+        call set_vertical(fu_set_layer_between_two(layer_btw_2_height, 1.0, 500.0), vertTmp)
         if(error)return
         arTmp(1) = 1.0
 
@@ -1012,7 +1013,8 @@ subroutine compute_emission_for_wb_dust(wb_dust_src, &
              & iMeteo, iDisp, ix, iy, ixMeteo, iyMeteo, iStat, iMode, iTmp, jTmp, iMeteoTmp
     real :: fTmp, timestep_sec, fCellTotal, ro_soil_dry, fSoilMassWater, u_star_min, &
           & uStarMerged, fSaltation, fFluxTotal, fRoughness, fSoilMassWaterThresh, fWindSpeed, &
-          & zx, zy, zz, fSum, du_star_owen, uStarGustEff !, uStarSaturated
+          & zx, zy, zz, fSum, du_star_owen, uStarGustEff, &
+          & soil_moisture, snow_depth, lai, fWindThreshold !, uStarSaturated
     real, dimension(:), pointer :: ptrXSzDisp, ptrYSzDisp, pSoilVolumeWater, &
                                  & fMinDArray, fMaxDArray, fPartDensity, pSandMassFract, &
                                  & pErodibleFraction, pU10m, pV10m, pU_star, pLAI, &
@@ -1044,7 +1046,7 @@ subroutine compute_emission_for_wb_dust(wb_dust_src, &
     !
     real, parameter :: fSandDensity = 2600, &
                      & fSandDensity_1_3 = fSandDensity ** 0.3333333, &
-                     & fSnowWEQDepthThreshold = 0.001, &     ! 1mm of snow => no dust
+                     & fSnowWEQDepthThreshold = 0.02, &     ! 2 cm of water-eq. snow -> no dust
                      & sqrt_density_air = sqrt(density_air_288K)
     logical, parameter :: ifWindGust = .true.
 
@@ -1222,7 +1224,6 @@ subroutine compute_emission_for_wb_dust(wb_dust_src, &
 fArDump = 0.0
 #endif
 
-
     !
     ! Now scan the domain, for each grid cell find the dust emission parameters
     ! and fill-in the emission map
@@ -1240,7 +1241,7 @@ subroutine compute_emission_for_wb_dust(wb_dust_src, &
           case(emis_sandblasting_v1_flag)
             if(pErodibleFraction(iDisp) < 1.0e-5)cycle
           case(simple_dust_flag)
-            if(pDustEmis0(iDisp) < 1.0e-18)cycle
+            if(pDustEmis0(iDisp) < 1.0e-14)cycle
           case default
         end select
 
@@ -1254,9 +1255,20 @@ subroutine compute_emission_for_wb_dust(wb_dust_src, &
 
 iArStatistics(2) = iArStatistics(2) + 1
 
+        if (pLandFr(iMeteo) < 0.95) then
+          call set_coastal_value(pSoilVolumeWater, pLandFr, iMeteo, nx_meteo, ny_meteo, soil_moisture, 0.3, .true.)
+          call set_coastal_value(pLAI, pLandFr, iMeteo, nx_meteo, ny_meteo, lai, 2.0, .false.)
+          call set_coastal_value(pSnowWEQDepth, pLandFr, iMeteo, nx_meteo, ny_meteo, snow_depth, 0.0, .false.)
+        else
+          soil_moisture = pSoilVolumeWater(iMeteo)
+          lai = pLAI(iMeteo)
+          snow_depth = pSnowWEQDepth(iMeteo)
+        end if
+
         !
         !   ...and if soil covered with snow
         !
+
         if(pSnowWEQDepth(iMeteo) > fSnowWEQDepthThreshold)cycle
 
 iArStatistics(3) = iArStatistics(3) + 1
@@ -1278,60 +1290,25 @@ subroutine compute_emission_for_wb_dust(wb_dust_src, &
           case default
             iArStatistics(5) = iArStatistics(5) + 1
         end select
-
-        !
-        ! Soil moisture can be tricky: for small islands and sea edges it can be set to zero. 
-        ! But it can also be zero for deserts. So, for partly sea-cells we will check the surrounding 
-        ! mainly land cells if they give any clue on actual moisture in the area
-        !
-        if(pLandFr(iMeteo) < 0.75)then
-          if(abs(pSoilVolumeWater(iMeteo)) < 1.0e-5)then
-if(ifTalking)call msg('Checking suspicious soil moisture:',ix,iy)
-            iCount = 0
-            fTmp = 0.0
-            do jTmp = max(iyMeteo-1,1), min(iyMeteo, ny_meteo)
-              do iTmp =  max(ixMeteo-1,1), min(ixMeteo, nx_meteo)
-                iMeteoTmp = iTmp + nx_meteo * (jTmp -1)
-                if(iMeteo == iMeteoTmp)cycle
-                if(pLandFr(iMeteoTmp) > 0.75)then
-                  fTmp = fTmp + pSoilVolumeWater(iMeteo)
-                  iCount = iCount + 1
-                endif
-              end do
-            end do
-            if(iCount > 0)then
-              fTmp = fTmp / real(iCount)
-              if(fTmp > 0.01)then
-                 pSoilVolumeWater(iMeteo) = fTmp  ! surrounding is wet: do not believe zero!
-if(ifTalking)call msg('Replace zero moisture at (' + fu_str(ix) + ',' + fu_str(iy) + ') with:', fTmp)
-              else
-if(ifTalking)call msg('Keep zero moisture')
-              endif
-            else
-if(ifTalking)call msg('Replace zero moisture at (' + fu_str(ix) + ',' + fu_str(iy) + ') with default 0.3')
-              pSoilVolumeWater(iMeteo) = 0.3  ! virtually anything: this is a kind-of isolated island
-            endif
-          endif ! soil moisture ~0
-        endif  ! significant part of water in cell 
 
         select case(wb_dust_src%emisMethod)
           case(simple_dust_flag)
 
             fWindSpeed = sqrt(pU10m(iMeteo)*pU10m(iMeteo)+pV10m(iMeteo)*pV10m(iMeteo))
-            fSaltation = max(1.4, fWindSpeed-4.1)**3 
-            fSaltation = fSaltation * max(0.0, ((1 - 5.*pSoilVolumeWater(iMeteo))))**2
+            fWindThreshold = 5 + 50*soil_moisture
+            fSaltation = max(1.4, fWindSpeed - fWindThreshold)**3
+
+            fSaltation = fSaltation * max(0.0, min(1.0, pLandFr(iMeteo)))
+            fSaltation = fSaltation * (max(0.0, fSnowWEQDepthThreshold - snow_depth)/fSnowWEQDepthThreshold)
 
             if(fSaltation <= 0) cycle
 
             ! pDustEmis0 refers to a precomputed scaling map that depends on the bare land fraction
             ! (i.e. source area mask), the surface rughness map and the sand fraction map.
-            ! Currently it equals scaling_factor * bare_land/roughness**0.5 * factor_based_on_land_features
-            ! 1/roughness**0.5 is numerically quite close to 1-log(roughness/some_roughness_scale),
-            ! but does not completely cut out the emission in rough areas (e.g. dust storms do occur
-            ! in the midwest USA). pDustEmis needed to be preprocessed to cut out mismatching map
-            ! pixels mainly along the coasts.
+            ! Currently it equals scaling_factor * bare_land_fraction * shape_function * ln(roughness/roughness_0)
+            ! The shape function represents alluvial deposits
 
-            fFluxTotal = 0.135 * pDustEmis0(iDisp) * fSaltation * timestep_sec * max(0.0, fLAI_critical-pLAI(iMeteo))
+            fFluxTotal = 0.178 * pDustEmis0(iDisp) * fSaltation * timestep_sec * max(0.0, fLAI_critical-lai)/fLAI_critical
 
             if(fSaltation > 0.0) iArStatistics(6) = iArStatistics(6) + 1
             if(fFluxTotal > 0.0) iArStatistics(7) = iArStatistics(7) + 1

source/toolbox.silja.mod.f90

@@ -90,6 +90,7 @@ MODULE toolbox
   public fu_get_default_mpi_buf_size
   public remapcon_1d
   public testremapcon_1d
+  public set_coastal_value
 
   public trim_precision
 
@@ -4915,7 +4916,64 @@ subroutine get_chunk(count_total, num_chunks, displs, sizes, if_zero_based)
 
   end subroutine get_chunk
 
-
+  !*************************************************************
+
+  subroutine set_coastal_value(field, pLandFr, iMeteo, nx_meteo, ny_meteo, &
+       & value_in_cell, default_value, if_divide)
+
+    real, dimension(:), pointer, intent(in) :: field, pLandFr
+    integer, intent(in) :: iMeteo, nx_meteo, ny_meteo
+    real, intent(in) :: default_value
+    logical, intent(in) :: if_divide
+    real, intent(out) :: value_in_cell
+    real :: fTmp
+    integer :: iCount, jTmp, iTmp, iyMeteo, ixMeteo, iMeteoTmp
+
+    ! Sets the value of a field in a coastal cell to equal the average of the
+    ! mostly land-containing neighbors. In case no mostly land-containing neighbor
+    ! is found, the default_value is used. The values can also be divided by the 
+    ! land area (controlled by if_divide).
+
+    if (pLandFr(iMeteo) < 0.85) then !if (field(iMeteo) < 1e-4) then
+      iCount = 0
+      fTmp = 0.0
+      do jTmp = max(iyMeteo-1,1), min(iyMeteo+1, ny_meteo)
+        do iTmp =  max(ixMeteo-1,1), min(ixMeteo+1, nx_meteo)
+          iMeteoTmp = iTmp + nx_meteo * (jTmp -1)
+          if (iMeteo == iMeteoTmp) cycle
+          if (pLandFr(iMeteoTmp) >= 0.85) then
+            if (if_divide) then
+              fTmp = fTmp + field(iMeteoTmp)/max(0.0, min(1.0, pLandFr(iMeteoTmp)))
+            else
+              fTmp = fTmp + field(iMeteoTmp)
+            end if
+            iCount = iCount + 1
+          end if
+        end do
+      end do
+      if (iCount > 0) then
+        fTmp = fTmp/iCount
+        if (fTmp > field(iMeteo)) then
+          value_in_cell = fTmp
+        else
+          if (if_divide) then
+            value_in_cell = field(iMeteo)/max(0.0, min(1.0, pLandFr(iMeteo)))
+          else
+            value_in_cell = field(iMeteo)
+          end if
+        endif
+      else
+        value_in_cell = default_value  ! virtually anything: this is a kind-of isolated island
+      endif
+    else
+      if (if_divide) then
+        value_in_cell = field(iMeteo)/max(0.0, min(1.0, pLandFr(iMeteo)))
+      else
+        value_in_cell = field(iMeteo)
+      end if
+    end if ! value ~0
+  end subroutine set_coastal_value
+
   !!!!!!***********************************************************************************
   !!!!!!
   !!!!!! A bunch of subroutines needed for IS4FIRES.