text:
!>
!! Initialization routine for multi-layer soil model
!!
!! <Describe the concepts of the procedures and algorithms used in the module.>
!! <Details of procedures are documented below with their definitions.>
!! <Include any applicable external references inline as module::procedure,>
!! <external_procedure(), or by using @see.>
!! <Don't forget references to literature.>
!!
!! @author Juergen Helmert, DWD
!! @author Ekaterina Machulskaya, DWD
!! @author Guenther Zaengl, DWD
!! @author Daniel Reinert, DWD
!!
!!
!! @par Revision History
!! Moved to separate module from mo_soil_ml.f90 by Daniel Reinert, DWED (2014-04-01)
!!
!! @par Copyright and License
!!
!! This code is subject to the DWD and MPI-M-Software-License-Agreement in
!! its most recent form.
!! Please see the file LICENSE in the root of the source tree for this code.
!! Where software is supplied by third parties, it is indicated in the
!! headers of the routines.
!!
MODULE mo_nwp_soil_init

  USE mo_kind,               ONLY: ireals=>wp, &
    &                              iintegers=>i4
  USE mo_run_config,         ONLY: msg_level
  USE mo_physical_constants, ONLY: t0_melt => tmelt,& ! absolute zero for temperature
    &                              rho_w => rhoh2o, & ! density of liquid water (kg/m^3)
    &                              lh_f  => alf   , & ! latent heat of fusion
    &                              g     => grav      ! acceleration due to gravity

  USE mo_lnd_nwp_config,     ONLY: lmulti_snow, lana_rho_snow, &
    &                              lmelt, lmelt_var,           &
    &                              max_toplaydepth, l2lay_rho_snow
  USE mo_phyparam_soil,      ONLY: cdw0, cdw1, ckw0, ckw1, cik2, &
    &                              cporv, cpwp, cadp, cfcap,     &
    &                              crock, crhoc, cala0, cala1,   &
    &                              csandf, cclayf, csnow_tmin,   &
   &                               crhosmax_tmin, crhosmax_ml

  IMPLICIT NONE

  PRIVATE

  PUBLIC  :: terra_multlay_init
  PUBLIC  :: get_wsnow

CONTAINS

  !>
  !! Performs initialization of soil model TERRA
  !!
  !! Performs initialization of soil model TERRA. Depending on the initialization mode
  !! chosen, a coldstart- , warmstart-  or warmstart with snow increments (IAU)-initialization
  !! is performed
  !!
  !!
  !! @par Revision History
  !! Initial Revision by Juergen Helmert, DWD (2011)
  !!
  SUBROUTINE terra_multlay_init (                &
                  init_mode,    & ! 1 = coldstart, 2 = warmstart, 3 = warmstart with snow increments (IAU)
                  ie,           & ! array dimensions
                  istartpar,    & ! start index for computations in the parallel program
                  iendpar,      & ! end index for computations in the parallel program
! 4. variables for the time discretization and related variables
!                  ke,           & ! nsubs0=1 for single tile, nsubs0=2 for multi-tile
!                  nsubs0ubs1,& ! nsubs1=1 for single tile, nsubs1=#tiles+1 for multi-tile
                  ke_soil, ke_snow , &
                  czmls            , & ! processing soil level structure
                  soiltyp_subs     , & ! type of the soil (keys 0-9)                     --
                  rootdp           , & ! depth of the roots                            ( m  )
                  plcov            , & ! fraction of surface covered by plants         ( -  )
                  t_snow_now       , & ! temperature of the snow-surface               (  K  )
                  t_snow_mult_now  , & ! temperature of the snow-surface               (  K  )
                  t_s_now          , & ! temperature of the ground surface             (  K  )
                  t_s_new          , & ! temperature of the ground surface             (  K  )
                  w_snow_now       , & ! water content of snow                         (m H2O)
                  h_snow           , & ! snow depth                                   (m H2O)
                  rho_snow_now     , & ! snow density                                  (kg/m**3)
                  rho_snow_mult_now, & ! snow density                                  (kg/m**3)
                  t_so_now         , & ! soil temperature (main level)                 (  K  )
                  t_so_new         , & ! soil temperature (main level)                 (  K  )
                  w_so_now         , & ! total water conent (ice + liquid water)       (m H20)
                  w_so_new         , & ! total water conent (ice + liquid water)       (m H20)
                  w_so_ice_now     , & ! ice content                                   (m H20)
                  w_so_ice_new     , & ! ice content                                   (m H20)
                  wliq_snow_now    , & ! liquid water content in the snow              (m H2O)
                  wtot_snow_now    , & ! total (liquid + solid) water content of snow  (m H2O)
                  dzh_snow_now       & ! layer thickness between half levels in snow   (  m  )
                                  )

!-------------------------------------------------------------------------------
! Declarations for ICON (USE statements for COSMO)
!-------------------------------------------------------------------------------

  IMPLICIT NONE

  INTEGER                 , INTENT(IN)  ::  &
                  init_mode    ! 1 = coldstart, 2 = warmstart, 3 = warmstart with snow increments (IAU)
  INTEGER (KIND=iintegers), INTENT(IN)  ::  &
                  ie,           & ! array dimensions
                  istartpar,    & ! start index for computations in the parallel program
                  iendpar,      & ! end index for computations in the parallel program
!                  ke,           &
!                 nsubs0ubs1 , &
                  ke_soil, ke_snow
  REAL    (KIND = ireals), DIMENSION(ke_soil+1), INTENT(IN) :: &
                  czmls           ! processing soil level structure

  INTEGER (KIND=iintegers),DIMENSION(ie), INTENT(IN) :: &
                  soiltyp_subs      ! type of the soil (keys 0-9)                     --
  REAL    (KIND = ireals), DIMENSION(ie), INTENT(IN) :: &
                  rootdp           ! depth of the roots                            ( m  )
  REAL    (KIND = ireals), DIMENSION(ie), INTENT(IN) :: &
                  plcov           ! plant coverage                            ( -  )
  REAL    (KIND = ireals), DIMENSION(ie), INTENT(INOUT) :: &
                  t_snow_now              ! temperature of the snow-surface (K)
  REAL    (KIND = ireals), DIMENSION(ie,0:ke_snow), INTENT(INOUT) :: &
                  t_snow_mult_now      ! temperature of the snow-surface               (  K  )
  REAL    (KIND = ireals), DIMENSION(ie), INTENT(INOUT) :: &
                  t_s_now              ! temperature of the ground surface             (  K  )
  REAL    (KIND = ireals), DIMENSION(ie), INTENT(OUT) :: &
                  t_s_new              ! temperature of the ground surface             (  K  )
  REAL    (KIND = ireals), DIMENSION(ie), INTENT(INOUT) :: &
                  w_snow_now       , & ! water content of snow                         (m H2O)
                  rho_snow_now         ! snow density                                  (kg/m**3)
  REAL    (KIND = ireals), DIMENSION(ie), INTENT(INOUT) :: &
                  h_snow               ! snow depth                                   (m H2O)
  REAL    (KIND = ireals), DIMENSION(ie,ke_snow), INTENT(INOUT) :: &
                  rho_snow_mult_now    ! snow density                                  (kg/m**3)
  REAL    (KIND = ireals), DIMENSION(ie,0:ke_soil+1), INTENT(INOUT) :: &
                  t_so_now             ! soil temperature (main level)                 (  K  )
  REAL    (KIND = ireals), DIMENSION(ie,0:ke_soil+1), INTENT(OUT) :: &
                  t_so_new             ! soil temperature (main level)                 (  K  )
  REAL    (KIND = ireals), DIMENSION(ie,ke_soil+1), INTENT(INOUT) :: &
                  w_so_now         , & ! total water conent (ice + liquid water)       (m H20)
                  w_so_ice_now         ! ice content                                   (m H20)
  REAL    (KIND = ireals), DIMENSION(ie,ke_soil+1), INTENT(OUT) :: &
                  w_so_new         , & ! total water conent (ice + liquid water)       (m H20)
                  w_so_ice_new         ! ice content                                   (m H20)
  REAL    (KIND = ireals), DIMENSION(ie,ke_snow), INTENT(INOUT) :: &
                  wliq_snow_now    , & ! liquid water content in the snow              (m H2O)
                  wtot_snow_now        ! total (liquid + solid) water content of snow  (m H2O)
  REAL    (KIND = ireals), DIMENSION(ie,ke_snow), INTENT(INOUT) :: &
                  dzh_snow_now         ! layer thickness between half levels in snow   (  m  )

!--------------------------------------------------------------------------------
! TERRA Declarations

! New declaration for ICON

!------------------------------------------------------------------------------
! Subroutine arguments: None
! --------------------

!
! Local parameters:
! ----------------

  REAL    (KIND=ireals   ), PARAMETER ::  &
    zepsi  = 1.0E-6_ireals ! security constant

  REAL    (KIND=ireals   ), PARAMETER ::  &
    T_ref_ice = 0.1_ireals , & !degC Soil ice parameterization
    T_star_ice= 0.01_ireals, & !degC according to K. Schaefer and Jafarov, E.,2016, doi:10.5194/bg-13-1991-2016
    b_clay= -0.3_ireals , &
    b_silt= -0.5_ireals , &
    b_sand= -0.9_ireals , &
    b_org = -1.0_ireals

! Local scalars:
! -------------

  INTEGER (KIND=iintegers) ::  &
    kso            , & ! loop index for soil moisture layers
    ksn,k          , & ! loop index for snow layers
    i,ic           , & ! loop index in x-direction
    icount_soil    , & ! "true" soil
    icount_rockice , & ! rock and ice points
    mstyp          , & ! soil type index
    istarts        , & ! start index for x-direction
    iends              ! end   index for x-direction

!  REAL    (KIND=ireals   ) ::  &
!    zdel_t_so          ! auxiliary variable

  REAL    (KIND=ireals   ) ::  &
    zpsi0=.01_ireals    ! air entry potential at water saturation (m)

  REAL    (KIND=ireals   ) ::  &
    zroota    (ie) ! root density profile parameter (1/m)

  REAL    (KIND=ireals   ) ::  &
    zaa, zh_snow, zzz, zd, zd1, zd2, zd3, zd4        ! utility variables

  REAL    (KIND = ireals) :: &
!
    zmls     (ke_soil+1)  , & ! depth of soil main level
    zzhls    (ke_soil+1)  , & ! depth of the half level soil layers in m
    zdzhs    (ke_soil+1)  , & ! layer thickness between half levels
    zdzms    (ke_soil+1)  , & ! distance between main levels

!
! External parameters
!
    zbwt     (ie)      , & ! root depth (with artificial minimum value)
    zsandf   (ie)      , & ! mean fraction of sand (weight percent)
    zclayf   (ie)      , & ! mean fraction of clay (weight percent)
    zsiltf   (ie)      , & ! mean fraction of silt (weight percent)
    zb_por   (ie)      , & ! pore size distribution index
    zpsis    (ie)      , & ! air entry potential (m)
    zw_m_org           , &  ! maximum of  liquid water content   organic
    zw_m_soil          , &  ! maximum of  liquid water content   mineral soil
    zw_m_up            , &  ! maximum of  liquid water content   at temp -3 degC
    zw_m_low           , &  ! maximum of  liquid water content   at temp -40 degC
    zw_m_org_low       , &  ! maximum of  liquid water content   organic at temp -40 degC
    zw_m_soil_low      , &  ! maximum of  liquid water content   mineral soil at temp -40 degC
    t_zw_up            , &  ! temp -3 degC
    t_zw_low           , &  ! temp -40 degC
    zw_m     (ie)          ! maximum of liquid water content  (m)

  INTEGER  (KIND=iintegers ) ::  &
    m_styp   (ie)      , & ! soil type
    zicount1 (ie)      , &
    zicount2 (ie)      , &
    soil_list(ie)      , & ! list of "true" soil points
                           ! i.e. no ice no rock
    rockice_list(ie)       ! list of rock and ice points
  REAL    (KIND=ireals   ) ::  & ! Soil and plant parameters
!
    zadp     (ie)      , & ! air dryness point
    zporv    (ie)      , & ! pore volume (fraction of volume)
    zedb     (ie)          ! utility variable

  REAL    (KIND=ireals   ) ::  &
    zw_snow_old(ie)      ,& !
    zrho_snow_old(ie)    ,&
    h_snow_fg(ie)        ,&
    h_snow_incr(ie)      ,&
    zhh_snow (ie,ke_snow),& ! depth of the half level snow layers
    zhm_snow (ie,ke_snow),& ! depth of snow main levels
    sum_weight(ie)      , &
    t_new  (ie,ke_snow) , &
    rho_new(ie,ke_snow) , &
    wl_new (ie,ke_snow) , &
    z_old  (ie,ke_snow) , &
    dz_old (ie,ke_snow) , &
    weight, wso_ice_tolerance, wso_ice_equil, rwso_ice_tolerance

  LOGICAL :: l_redist(ie)

!- End of header
!==============================================================================

!------------------------------------------------------------------------------
! Begin Subroutine terra_multlay
!------------------------------------------------------------------------------
!==============================================================================
!  Computation of the diagnostic part I of the soil parameterization scheme
!  In this part, evaporation from the surface and transpiration is calculated.
!  A multi-layer soil water distribution is calculated by simple bulk
!  parameterisation or a Penman-Monteith-type version of evaporation
!  and transpiration which can be used alternatively.
!------------------------------------------------------------------------------

!------------------------------------------------------------------------------
! Section I.1: Initializations
!------------------------------------------------------------------------------

! Horizontal domain for computation
  istarts = istartpar
  iends   = iendpar

! grids for temperature and water content

  zzhls(1) = 2._ireals*czmls(1)   !depth of first half level
  zdzhs(1) = zzhls(1)      !layer thickness betw. half levels of uppermost layer
  zmls(1)  = czmls(1)      !depth of 1st main level
  zdzms(1) = czmls(1)      !layer thickness between soil surface and main level
                           ! of uppermost layer

  DO kso = 2,ke_soil+1
    zzhls(kso)  = zzhls(kso-1) + 2._ireals*(czmls(kso) -zzhls(kso-1))
    zdzhs(kso) = zzhls(kso) - zzhls(kso-1) ! layer thickness betw. half levels
    zmls(kso)  = czmls(kso)                ! depth of main levels
    zdzms(kso) = zmls(kso) - zmls(kso-1)   ! layer thickness betw. main levels
  ENDDO

! Prepare basic surface properties (for land-points only)
  icount_soil=0
  icount_rockice=0
  DO i = istarts, iends
    mstyp       = soiltyp_subs(i)        ! soil type
    m_styp(i) = mstyp                     ! array for soil type
    IF (mstyp >= 3) THEN
      icount_soil=icount_soil+1
      soil_list(icount_soil)=i
    ELSE
      icount_rockice=icount_rockice+1
      rockice_list(icount_rockice)=i
    END IF
    zporv(i)  = cporv(mstyp)              ! pore volume
    zadp (i)  = cadp (mstyp)              ! air dryness point
    zbwt(i)   = MAX(0.001_ireals,rootdp(i))! Artificial minimum value
                                           ! for root depth
    zroota(i) = 3._ireals/zbwt(i)       ! root density profile parameter (1/m)
                                        ! zroota=0. creates the original TERRA_LM
                                        ! version with constant root density
                                        ! for root depth

    ! Arrays for soil water freezing/melting
    zsandf(i)   = csandf(mstyp)
    zclayf(i)   = cclayf(mstyp)
    zsiltf(i)   = 100._ireals -csandf(mstyp)-cclayf(mstyp) ! Residuum of sand and clay
    zpsis(i)    = -zpsi0 * 10._ireals**(1.88_ireals-0.013_ireals*zsandf(i))
    zb_por(i)   = 2.91_ireals + .159_ireals*zclayf(i)
    zedb(i)     = 1._ireals/zb_por(i)

! print*,i,zporv(i),zclayf(i),zpsis(i),zb_por(i)

  ENDDO

  ! Further parameters for soil water freezing/melting
  t_zw_up  = 270.15_ireals ! temp -3 degC
  t_zw_low = 233.15_ireals ! temp -40 degC
  zd = LOG((T_ref_ice-(t_zw_low-t0_melt))/T_star_ice)
  zd1 = EXP(b_sand*zd)
  zd2 = EXP(b_clay*zd)
  zd3 = EXP(b_silt*zd)
  zd4 = EXP(b_org*zd)

  ! Tolerances for the allowed deviation of w_so_ice from its equilibrium value
  wso_ice_tolerance  = 1.05_ireals ! 5% deviation allowed
  rwso_ice_tolerance = 1._ireals/wso_ice_tolerance

! For ntstep=0 : Some preparations
! ================================

  w_so_new(:,:) = w_so_now(:,:)

! SBr, CHa: DO NOT set minimal values to w_so if the threshold is passed
!! Provide for a soil moisture 1 % above air dryness point, reset soil
!! moisture to zero in case of ice and rock
!  DO kso   = 1,ke_soil+1
!    DO ic = 1, icount_soil
!      i=soil_list(ic)
!      w_so_now (i,kso) = MAX(w_so_now(i,kso),                     &
!                         1.01_ireals*zadp(i)*zdzhs(kso) )
!      w_so_new (i,kso) = MAX(w_so_new(i,kso),                     &
!                         1.01_ireals*zadp(i)*zdzhs(kso) )
!    END DO
!    DO ic = 1, icount_rockice
!      i=rockice_list(ic)
!      w_so_now(i,kso) = 0.0_ireals
!      w_so_new(i,kso) = 0.0_ireals
!    END DO
!  END DO

! adjust temperature profile in lower soil layers, if temperature of first soil
! layer was reduced due to the detection of snow (e.g. in the analysis)
! loop over grid points
  DO i = istarts, iends
    IF (w_snow_now(i) <= zepsi) THEN
      ! spurious snow is removed
      t_snow_now(i)=t_so_now(i,0)
!jh   t_snow_new(i)=t_so_now(i,0)
      w_snow_now(i)= 0.0_ireals
      IF(lmulti_snow) THEN
        t_snow_mult_now(i,0) = t_so_now(i,0)
        DO ksn = 1, ke_snow
          t_snow_mult_now(i,ksn) = t_so_now(i,0)
          wliq_snow_now(i,ksn) = 0.0_ireals
          wtot_snow_now(i,ksn) = 0.0_ireals
          rho_snow_mult_now(i,ksn) = 0.0_ireals
          dzh_snow_now(i,ksn) = 0.0_ireals
        END DO
      END IF
! GZ: this soil temperature adjustment tends to cause more damage than benefit in practice
! because if there are traces of snow over warm ground, soil heat fluxes can become so large
! that the scheme becomes numerically unstable if the soil top temperature is artificially
! adjusted to freezing temperature
    ELSE
! !           adjust soil temperatures in presence of snow
!             zdel_t_so = MAX (0._ireals,t_so_now(i,0)-t0_melt)
!             t_so_now(i,0)=MIN( t0_melt,t_so_now(i,0) )
!             t_so_new(i,0)=MIN( t0_melt,t_so_now(i,0) )
!             DO kso=1,ke_soil
!               IF ( t_so_now(i,kso) > t0_melt) THEN
!                  t_so_now(i,kso) = MAX( t0_melt,                  &
!                                    t_so_now(i,kso) -zdel_t_so    &
!                                    *(zmls(ke_soil+1)-zmls(kso))   &
!                                    /(zmls(ke_soil+1)-zmls( 1 )) )
!                  t_so_new(i,kso) = MAX( t0_melt,                  &
!                                    t_so_now(i,kso) -zdel_t_so    &
!                                    *(zmls(ke_soil+1)-zmls(kso))   &
!                                    /(zmls(ke_soil+1)-zmls( 1 )) )
!               ENDIF
!               IF ( t_so_now(i,kso) <= t0_melt) EXIT
!             ENDDO
      t_snow_now(i) = MIN( t0_melt, t_snow_now(i))
      IF(lmulti_snow) THEN
        DO ksn = 0, ke_snow
          t_snow_mult_now(i,ksn) = MIN( t0_melt, t_snow_mult_now(i,ksn))
        END DO
      END IF
    ENDIF

    t_s_now(i) = t_so_now(i,0)
    t_s_new(i) = t_so_now(i,0)

!   Set level 1 to level 0 for t_so for every landpoint
    t_so_now(i,1) = t_so_now(i,0)
    t_so_new(i,1) = t_so_now(i,0)
    t_so_new(i,0) = t_so_now(i,0)

  END DO

  ! Note that a few parts need to be shifted to the .NOT. is_coldstart branch.
  ! I.e. those which are dealing with the contributions from the analysis.
  ! Currently, the multi layer soil model starts from the FG.
  IF ( init_mode == 1 ) THEN

    ! Initialization of
    ! - snow density (if necessary)
    ! - multi layer snow variables (if necessary)
    !   * wtot_snow, dzh_snow, wliq_snow
    ! - soil ice water content w_so_ice
    ! - snow height h_snow

!   Initialization of snow density, if necessary
!   --------------------------------------------
    IF (.NOT. lana_rho_snow) THEN
       IF(lmulti_snow) THEN
         ksn = 1
         DO i = istarts, iends
           rho_snow_mult_now(i,ksn) = 250.0_ireals    ! average initial density
           t_snow_mult_now  (i,ksn) = t_snow_now(i)
           wliq_snow_now    (i,ksn) = 0.0_ireals
           ! Limit top layer to max_toplaydepth
           dzh_snow_now     (i,ksn) = MIN(max_toplaydepth, w_snow_now(i)/REAL(ke_snow,ireals) &
                  &                            *rho_w/rho_snow_mult_now(i,ksn))
           wtot_snow_now    (i,ksn) = dzh_snow_now(i,ksn)*rho_snow_mult_now(i,ksn)/rho_w
           t_snow_mult_now  (i,0  ) = t_snow_now(i)
         ENDDO
         k = MIN(2, ke_snow-1)
         DO ksn = 2, ke_snow
           DO i = istarts, iends
             rho_snow_mult_now(i,ksn) = 250.0_ireals    ! average initial density
             t_snow_mult_now  (i,ksn) = t_snow_now(i)
             wliq_snow_now    (i,ksn) = 0.0_ireals
             ! distribute the remaining snow equally among the layers
             IF (k == ksn) THEN
               dzh_snow_now   (i,ksn) = MIN(8._ireals*max_toplaydepth, (w_snow_now(i)-wtot_snow_now(i,1)) &
                    &                            /REAL(ke_snow-1,ireals)*rho_w/rho_snow_mult_now(i,ksn) )
             ELSE
               dzh_snow_now   (i,ksn) = (w_snow_now(i)-wtot_snow_now(i,k))/REAL(ke_snow-k,ireals) &
                    &                            *rho_w/rho_snow_mult_now(i,ksn)
             ENDIF
             wtot_snow_now    (i,ksn) = dzh_snow_now(i,ksn)*rho_snow_mult_now(i,ksn)/rho_w
           ENDDO
         ENDDO
       ELSE
         rho_snow_now(istarts:iends) = 250.0_ireals    ! average initial density
         IF (l2lay_rho_snow) rho_snow_mult_now(istarts:iends,:) = 250.0_ireals
       ENDIF
    ELSE
      IF(lmulti_snow) THEN
!The sum of wtot_snow, i.e. the "old" total water equivalent depth
        zw_snow_old  (:) = 0.0_ireals
        zrho_snow_old(:) = 0.0_ireals
        zicount1(:)      = 0
        zicount2(:)      = 0
        sum_weight(:)    = 0.0_ireals
        DO ksn = 1, ke_snow
          DO i = istarts, iends
            zw_snow_old(i) = zw_snow_old(i) + wtot_snow_now(i,ksn)
            zrho_snow_old(i) = zrho_snow_old(i) + dzh_snow_now(i,ksn)
            IF (wtot_snow_now(i,ksn) > zepsi) zicount1(i) = zicount1(i)+1
            IF (dzh_snow_now(i,ksn) > zepsi)  zicount2(i) = zicount2(i)+1
          END DO
        END DO
!The "old" snow density
        DO i = istarts, iends
          zrho_snow_old(i) = zw_snow_old(i) / MAX(zrho_snow_old(i),1.0E-09_ireals) *rho_w
        END DO
        k = MIN(2, ke_snow-1)
        DO ksn = 1, ke_snow
          DO i = istarts, iends
            IF(zicount1(i).EQ.ke_snow .AND. zicount2(i).EQ.ke_snow) THEN
              rho_snow_mult_now(i,ksn) = wtot_snow_now(i,ksn)/dzh_snow_now(i,ksn)*rho_w
              wtot_snow_now(i,ksn) = wtot_snow_now(i,ksn)*w_snow_now(i)/zw_snow_old(i)
              dzh_snow_now (i,ksn) = dzh_snow_now (i,ksn)*w_snow_now(i)/zw_snow_old(i)
              wliq_snow_now(i,ksn) = wliq_snow_now(i,ksn)*w_snow_now(i)/zw_snow_old(i)
            ELSE
              IF(rho_snow_now(i) .EQ. 0._ireals) rho_snow_now(i) = 250._ireals
              IF (ksn == 1) THEN ! Limit top layer to max_toplaydepth
                IF (t_snow_now(i) > 268.15_ireals) THEN
                  rho_snow_mult_now(i,ksn) = MIN(400._ireals,rho_snow_now(i))
                ELSE ! at very low temperatures, the near-surface snow is usually not subject to strong densification
                  rho_snow_mult_now(i,ksn) = MIN(rho_snow_now(i),400._ireals+6._ireals* &
                    (MAX(228.15_ireals,t_snow_now(i))-268.15_ireals)                    )
                ENDIF
                dzh_snow_now(i,ksn) = MIN(max_toplaydepth, w_snow_now(i)/REAL(ke_snow,ireals) &
                  &                  *rho_w/rho_snow_mult_now(i,ksn))
              ELSE ! distribute the remaining snow equally among the layers
                rho_snow_mult_now(i,ksn) = rho_snow_now(i)
                IF (ksn == k) THEN
                  dzh_snow_now(i,ksn) = MIN(8._ireals*max_toplaydepth, (w_snow_now(i)-wtot_snow_now(i,1)) &
                      &                 /REAL(ke_snow-1,ireals)*rho_w/rho_snow_mult_now(i,ksn) )
                ELSE
                  dzh_snow_now(i,ksn) = (w_snow_now(i)-sum_weight(i))/REAL(ke_snow-k,ireals) &
                      &                 *rho_w/rho_snow_mult_now(i,ksn)
                ENDIF
              ENDIF
              wtot_snow_now    (i,ksn) = dzh_snow_now(i,ksn)*rho_snow_mult_now(i,ksn)/rho_w
              sum_weight(i) = sum_weight(i) + wtot_snow_now(i,ksn)
              wliq_snow_now    (i,ksn) = 0.0_ireals
            END IF
          END DO
        END DO
        DO i = istarts, iends
          t_snow_mult_now  (i,0  ) = t_snow_now(i)
        END DO
      ELSE
        DO i = istarts, iends
          IF(rho_snow_now(i) .EQ. 0._ireals) rho_snow_now(i) = 250._ireals
          IF (l2lay_rho_snow) THEN
            IF(rho_snow_mult_now(i,1) == 0._ireals) rho_snow_mult_now(i,1) = rho_snow_now(i)
            IF(rho_snow_mult_now(i,2) == 0._ireals) rho_snow_mult_now(i,2) = rho_snow_now(i)
          ENDIF
        END DO
      END IF
    ENDIF

!   Initialization of soil ice content
!   ----------------------------------
!   Generally the combination lmelt=.true., lmelt_var=.true. has to be used.
!   lmelt_var=.false. (soil water freezing at 273.15K) should be used only
!   for special investigations because at model start a partial frozen layer
!   is only considered if the soil ice water content is read in from
!   a pre run.
!   ----------------------------------
    IF (lmelt .AND. .NOT. lmelt_var) THEN
      DO kso   = 1,ke_soil+1
        DO i = istarts, iends
          w_so_ice_now(i,kso) = 0.0_ireals
          w_so_ice_new(i,kso) = 0.0_ireals

          IF (t_so_now(i,kso) < t0_melt) THEN
            w_so_ice_now(i,kso) = w_so_now(i,kso)
            w_so_ice_new(i,kso) = w_so_now(i,kso)
          END IF ! t_so(kso) < t0_melt
        END DO
      END DO
    END IF           ! lmelt .AND. .NOT. lmelt_var
    IF(lmelt .AND. lmelt_var) THEN
      DO kso   = 1,ke_soil+1
        DO i = istarts, iends
          IF (t_so_now(i,kso) < (t0_melt-zepsi)) THEN

            zaa    = g*zpsis(i)/lh_f

            ! J. Helmert: Soil ice parameterization according to K. Schaefer and Jafarov, E.,2016,
            !                                                    doi:10.5194/bg-13-1991-2016

            ! Liq. water content at -3 degC
            zw_m_up =  zporv(i)*zdzhs(kso)*EXP(-zedb(i)*LOG((t_zw_up - t0_melt)/(t_zw_up*zaa)) )

            ! Determine liq. water content at -40 degC
            zw_m_soil = 0.01_ireals*(zsandf(i)*zd1 + zclayf(i)*zd2 + zsiltf(i)*zd3)

            ! Scale soil ice content with organic soil horizon.
            ! should decrease the root zone liquid water content of frozen soil for low temperatures significantly!
            IF(zmls(kso) < rootdp(i)) THEN
              zzz = plcov(i)*(rootdp(i)-zmls(kso))/rootdp(i)
              zw_m_low = zporv(i)*zdzhs(kso)*(zzz*zd4 + (1._ireals-zzz)*zw_m_soil)
            ELSE
              zw_m_low = zporv(i)*zdzhs(kso)*zw_m_soil
            END IF

            IF (t_so_now(i,kso) < t_zw_low) THEN
              zw_m(i) = zw_m_low
            ELSE IF (t_so_now(i,kso) < t_zw_up) THEN ! Logarithmic Interpolation between -3 degC and -40 degC
              zw_m(i) = zw_m_low*EXP((t_so_now(i,kso) - t_zw_low)*(LOG(zw_m_up) - LOG(zw_m_low))/(t_zw_up-t_zw_low))
            ELSE
              zw_m(i) = zporv(i)*zdzhs(kso)*EXP(-zedb(i)*LOG((t_so_now(i,kso)-t0_melt)/(t_so_now(i,kso)*zaa)) )
            END IF

            w_so_ice_now(i,kso) = MAX (0.0_ireals,w_so_now(i,kso) - zw_m(i))
            w_so_ice_new(i,kso) = MAX (0.0_ireals,w_so_now(i,kso) - zw_m(i))
          ELSE ! ensure that w_so_ice is zero
            w_so_ice_now(i,kso) = 0.0_ireals
            w_so_ice_new(i,kso) = 0.0_ireals
          END IF
        END DO
      END DO
    END IF           ! lmelt .AND. lmelt_var

!   Initialization of the local array of the grid in snow
!   -----------------------------------------------------
    IF(lmulti_snow) THEN
      sum_weight(:)    = 0.0_ireals
      DO i = istarts, iends
        h_snow(i) = 0.0_ireals
        DO ksn = 1,ke_snow
          h_snow(i) = h_snow(i) + dzh_snow_now(i,ksn)
        END DO
      END DO
      DO ksn = 1,ke_snow
        DO i = istarts, iends
          IF(h_snow(i) > 0.0_ireals) THEN
            t_snow_mult_now(i,ksn) = t_snow_now(i) + &
              (t_s_now(i)-t_snow_now(i))*(sum_weight(i)+dzh_snow_now(i,ksn)*0.5)/h_snow(i)
            sum_weight(i) = sum_weight(i) + dzh_snow_now(i,ksn)
          END IF
        END DO
      END DO
    ELSE
      DO i = istarts, iends
        h_snow(i) = w_snow_now(i)/rho_snow_now(i)*rho_w
      END DO
    END IF

  ELSE  ! init_mode > 1, i.e. warmstart (i.e. assimilation cycle and forecast)
        ! in this case, w_snow and rho_snow are read from the first guess, and h_snow is updated by the snow analysis;
        ! in the IAU mode, the snow depth increment (h_snow_incr) is directly provided by the snow analysis

    IF ( init_mode == 2 ) THEN  ! snow depth increment needs to be computed
      DO i = istarts, iends
        h_snow_fg(i)   = w_snow_now(i)/rho_snow_now(i)*rho_w
        ! Self-consistency fix needed because of possible GRIB truncation errors
        IF (lmulti_snow .AND. h_snow_fg(i) > zepsi) THEN
          dzh_snow_now(i,2:ke_snow) = MAX(dzh_snow_now(i,1),dzh_snow_now(i,2:ke_snow))
          zh_snow       = SUM(dzh_snow_now(i,1:ke_snow))
          w_snow_now(i) = w_snow_now(i) * zh_snow/h_snow_fg(i)
          h_snow_fg(i)  = zh_snow
        ELSE IF (lmulti_snow) THEN
          dzh_snow_now(i,:) = 0._ireals
        ENDIF
        h_snow_incr(i) = h_snow(i) - h_snow_fg(i)
      ENDDO
    ELSE IF ( init_mode == 3 ) THEN  ! snow depth increment is provided from snow analysis
      IF (lmulti_snow) THEN
        DO i = istarts, iends
          ! Self-consistency fix needed because of possible GRIB truncation errors
          IF (dzh_snow_now(i,1) > zepsi) THEN
            dzh_snow_now(i,2:ke_snow) = MAX(dzh_snow_now(i,1),dzh_snow_now(i,2:ke_snow))
            h_snow_fg(i) = SUM(dzh_snow_now(i,1:ke_snow))
          ELSE
            dzh_snow_now(i,:) = 0._ireals
            h_snow_fg(i)      = 0._ireals
          ENDIF
          h_snow_incr(i) = h_snow(i) - h_snow_fg(i)
        ENDDO
      ELSE
        DO i = istarts, iends
          h_snow_fg(i) = h_snow(i)
          h_snow_incr(i) = 0._ireals
        ENDDO
      ENDIF
    ENDIF

    ! Provisional fix to get rid of isolated extremely cold points in parallel routine (2014-10-16)
    ! This is applied only to the lowest prognostic soil layer, i.e. ke_soil
    DO i = istarts, iends
      IF (MIN(t_so_now(i,ke_soil+1),t_so_now(i,ke_soil-1)) - t_so_now(i,ke_soil) > 5._ireals) THEN
        t_so_now(i,ke_soil) = MIN(t_so_now(i,ke_soil+1),t_so_now(i,ke_soil-1)) - 5._ireals
      ENDIF
      IF (t_so_now(i,ke_soil) - MAX(t_so_now(i,ke_soil+1),t_so_now(i,ke_soil-1)) > 5._ireals) THEN
        t_so_now(i,ke_soil) = MAX(t_so_now(i,ke_soil+1),t_so_now(i,ke_soil-1)) + 5._ireals
      ENDIF
    ENDDO

    ! Ensure that w_so_ice stays within 5% of its equilibrium value
    ! In addition, it must not exceed w_so
    IF(lmelt .AND. lmelt_var) THEN
      DO kso   = 1,ke_soil+1
        DO i = istarts, iends
          IF (t_so_now(i,kso) < (t0_melt-zepsi)) THEN

            zaa    = g*zpsis(i)/lh_f

            ! J. Helmert: Soil ice parameterization according to K. Schaefer and Jafarov, E.,2016,
            !                                                    doi:10.5194/bg-13-1991-2016

            ! Liq. water content at -3 degC
            zw_m_up =  zporv(i)*zdzhs(kso)*EXP(-zedb(i)*LOG((t_zw_up - t0_melt)/(t_zw_up*zaa)) )

            ! Determine liq. water content at -40 degC
            zw_m_soil = 0.01_ireals*(zsandf(i)*zd1 + zclayf(i)*zd2 + zsiltf(i)*zd3)

            ! Scale soil ice content with organic soil horizon.
            ! should decrease the root zone liquid water content of frozen soil for low temperatures significantly!
            IF(zmls(kso) < rootdp(i)) THEN
              zzz = plcov(i)*(rootdp(i)-zmls(kso))/rootdp(i)
              zw_m_low = zporv(i)*zdzhs(kso)*(zzz*zd4 + (1._ireals-zzz)*zw_m_soil)
            ELSE
              zw_m_low = zporv(i)*zdzhs(kso)*zw_m_soil
            END IF

            IF (t_so_now(i,kso) < t_zw_low) THEN
              zw_m(i) = zw_m_low
            ELSE IF (t_so_now(i,kso) < t_zw_up) THEN ! Logarithmic Interpolation between -3 degC and -40 degC
              zw_m(i) = zw_m_low*EXP((t_so_now(i,kso) - t_zw_low)*(LOG(zw_m_up) - LOG(zw_m_low))/(t_zw_up-t_zw_low))
            ELSE
              zw_m(i) = zporv(i)*zdzhs(kso)*EXP(-zedb(i)*LOG((t_so_now(i,kso)-t0_melt)/(t_so_now(i,kso)*zaa)) )
            END IF

            wso_ice_equil = MAX (0.0_ireals,w_so_now(i,kso) - zw_m(i))
            w_so_ice_now(i,kso) = MIN(w_so_ice_now(i,kso),  wso_ice_tolerance*wso_ice_equil, w_so_now(i,kso))
            w_so_ice_now(i,kso) = MAX(w_so_ice_now(i,kso), rwso_ice_tolerance*wso_ice_equil)
          ELSE ! ensure that w_so_ice is zero
            w_so_ice_now(i,kso) = 0.0_ireals
          END IF
        END DO
      END DO
    END IF

    ! For the single-layer snow case, the prognostic variables now need to be updated. This requires appropriate
    ! assumptions about the new snow density

    IF(lmulti_snow) THEN ! update of prognostic variables of multi-layer snow model
      l_redist(:) = .FALSE.
      DO i = istarts, iends
        IF (h_snow(i) <= zepsi) THEN  ! snow has disappeared or was already absent
          rho_snow_now(i) = 250._ireals
          w_snow_now(i)   = 0._ireals

          dzh_snow_now (i,1:ke_snow) = 0._ireals
          wtot_snow_now(i,1:ke_snow) = 0._ireals
          rho_snow_mult_now(i,1:ke_snow) = 0._ireals
          wliq_snow_now(i,1:ke_snow) = 0._ireals
          t_snow_mult_now(i,1:ke_snow) = t_snow_now(i)
        ELSE IF (h_snow_fg(i) <= zepsi .AND. h_snow(i) > zepsi) THEN

          ! snow analysis has either reestablished a snow cover that had been erroneously melted away by the model
          ! or generated a new snow cover that was missed by the model due to a lack of precipitation (or wrong phase of precip)
          ! ** snow density should be distinguished between new snow and reestablished old snow **
          rho_snow_now(i) = 250._ireals ! needs to be improved
          w_snow_now(i)   = h_snow(i)*rho_snow_now(i)/rho_w
          t_snow_now(i)   = MIN(t_snow_now(i),t0_melt)

          rho_snow_mult_now(i,1:ke_snow) = rho_snow_now(i)
          dzh_snow_now (i,1:ke_snow) = h_snow(i)/REAL(ke_snow,ireals)
          wtot_snow_now(i,1:ke_snow) = w_snow_now(i)/REAL(ke_snow,ireals)
          wliq_snow_now(i,1:ke_snow) = 0._ireals
          t_snow_mult_now(i,1:ke_snow) = t_snow_now(i)
          IF (h_snow(i) > REAL(ke_snow,ireals)*max_toplaydepth) l_redist(i) = .TRUE.
        ELSE IF (h_snow_fg(i) > zepsi .AND. h_snow_incr(i) < 0._ireals) THEN

          ! a snow cover is present but reduced by the snow analysis increment
          ! in this case, we retain the snow densities but adjust the water equivalents
          l_redist(i) = .TRUE.
          dz_old(i,1) = MAX(dzh_snow_now (i,1),zepsi/REAL(ke_snow,ireals))
          dzh_snow_now (i,1) = dzh_snow_now (i,1) + h_snow_incr(i)
          wtot_snow_now(i,1) = wtot_snow_now(i,1)*dzh_snow_now (i,1)/dz_old(i,1)
        ELSE
          ! a snow cover is present and increased by the snow analysis
          ! ** snow density should be distinguished between new snow and reestablished old snow **
          ! IF (w_snow_incr(i) > 0._ireals) THEN ! assume temperature-dependent fresh snow density - maybe air temp should be used here?!
          !   rho_new(i,1) = MAX(50._ireals,125._ireals-10._ireals*(t0_melt-t_snow_mult_now(i,1)))
          ! ELSE
          !   rho_new(i,1) = rho_snow_mult_now(i,1)
          ! ENDIF
          l_redist(i) = .TRUE.
          rho_new(i,1) = rho_snow_mult_now(i,1)
          wtot_snow_now(i,1) = wtot_snow_now(i,1) + h_snow_incr(i)*rho_new(i,1)/rho_w
          dzh_snow_now (i,1) = dzh_snow_now (i,1) + h_snow_incr(i)
          rho_snow_mult_now(i,1) = rho_w*wtot_snow_now(i,1)/dzh_snow_now (i,1)
        ENDIF
      ENDDO

      ! Redistribute snow cover among layers where necessary
      sum_weight(:) = 0._ireals
      DO ksn = ke_snow,1,-1
        DO i = istarts, iends
          IF (l_redist(i)) THEN
            dz_old(i,ksn) = dzh_snow_now(i,ksn)
            z_old(i,ksn) = -sum_weight(i) - dzh_snow_now(i,ksn)/2._ireals
            sum_weight(i) = sum_weight(i) + dzh_snow_now(i,ksn)
          ENDIF
        END DO
      END DO

      k = MIN(2, ke_snow-1)
      DO ksn = 1,ke_snow
        DO i = istarts, iends
          IF (l_redist(i)) THEN

            IF (ksn == 1) THEN ! Limit top layer to max_toplaydepth
              zhh_snow(i,ksn) = -MAX( h_snow(i)-max_toplaydepth, h_snow(i)/ke_snow*(ke_snow-ksn) )
              zhm_snow(i,ksn) = (-h_snow(i) + zhh_snow(i,ksn))/2._ireals
              dzh_snow_now(i,ksn) = zhh_snow(i,ksn) + h_snow(i)    !layer thickness betw. half levels of uppermost snow layer
            ELSE IF (ksn == 2 .AND. ke_snow > 2) THEN ! Limit second layer to 8*max_toplaydepth
              zhh_snow(i,ksn) = MIN( 8._ireals*max_toplaydepth+zhh_snow(i,1), zhh_snow(i,1)/(ke_snow-1)*(ke_snow-ksn) )
            ELSE ! distribute the remaining snow equally among the layers
              zhh_snow(i,ksn) = zhh_snow(i,k)/(ke_snow-k)*(ke_snow-ksn)
            ENDIF

            t_new  (i,ksn) = 0.0_ireals
            rho_new(i,ksn) = 0.0_ireals
            wl_new (i,ksn) = 0.0_ireals

            IF(dz_old(i,ksn) > 0._ireals .AND. rho_snow_mult_now(i,ksn) > 0._ireals) THEN
              wliq_snow_now(i,ksn) = wliq_snow_now(i,ksn)/dz_old(i,ksn)
            END IF

          ENDIF
        END DO
      END DO

      DO ksn = 2, ke_snow
        DO i = istarts, iends
          IF(l_redist(i)) THEN
            zhm_snow(i,ksn) = (zhh_snow(i,ksn) + zhh_snow(i,ksn-1))/2._ireals
            dzh_snow_now(i,ksn) = zhh_snow(i,ksn) - zhh_snow(i,ksn-1) ! layer thickness betw. half levels
          END IF
        END DO
      END DO

      ! mass-weighted recomputation of temperatures and densities
      DO ksn = ke_snow,1,-1
        DO k = ke_snow,1,-1
          DO i = istarts, iends
            IF(l_redist(i)) THEN

              weight = MAX(MIN(z_old(i,k)+dz_old(i,k)/2._ireals,zhm_snow(i,ksn) + dzh_snow_now(i,ksn)/2._ireals)-   &
                       MAX(z_old(i,k)-dz_old(i,k)/2._ireals, zhm_snow(i,ksn)-dzh_snow_now(i,ksn)/2._ireals),0._ireals) &
                       &/dzh_snow_now(i,ksn)

              t_new  (i,ksn) = t_new  (i,ksn) + t_snow_mult_now  (i,k)*weight
              rho_new(i,ksn) = rho_new(i,ksn) + rho_snow_mult_now(i,k)*weight
              wl_new (i,ksn) = wl_new (i,ksn) + wliq_snow_now    (i,k)*weight
            END IF
          END DO
        END DO
      END DO

      DO ksn = ke_snow,1,-1
        DO i = istarts, iends
          IF(l_redist(i)) THEN
            t_snow_mult_now  (i,ksn) = t_new  (i,ksn)
            rho_snow_mult_now(i,ksn) = rho_new(i,ksn)
            wtot_snow_now    (i,ksn) = rho_new(i,ksn)*dzh_snow_now(i,ksn)/rho_w
            wliq_snow_now    (i,ksn) = wl_new (i,ksn)*dzh_snow_now(i,ksn)
          END IF
        END DO
      END DO

      ! Re-diagnose integrated/averaged snow-fields, in order to
      ! have a consistent state.
      !
      h_snow      (istarts:iends) = 0.0_ireals
      w_snow_now  (istarts:iends) = 0.0_ireals
      rho_snow_now(istarts:iends) = 0.0_ireals
      t_snow_now  (istarts:iends) = 0.0_ireals
      ! Re-diagnose h_snow, w_snow, rho_snow, t_snow
      DO ksn = ke_snow,1,-1
        DO i = istarts, iends
          h_snow(i)    = h_snow(i)     + dzh_snow_now(i,ksn)
          w_snow_now(i)= w_snow_now(i) + wtot_snow_now(i,ksn)
        END DO
      ENDDO

      DO ksn = ke_snow,1,-1
        DO i = istarts, iends
          IF (h_snow(i) > zepsi) THEN
            rho_snow_now(i) = rho_snow_now(i) + rho_snow_mult_now(i,ksn) * dzh_snow_now(i,ksn)/h_snow(i)
          ELSE
            rho_snow_now(i) = 250._ireals
          ENDIF
        END DO
      ENDDO

      DO i = istarts, iends
        t_snow_now(i) = t_snow_mult_now(i,1)
      ENDDO

    ELSE  ! single snow layer

      ! diagnose w_snow and rho_snow from analyzed h_snow
      ! h_snow is taken from analysis, while up to this point
      ! rho_snow_now and w_snow_now contain the first guess.
      CALL get_wsnow(h_snow,         &  ! in
        &            rho_snow_now,   &  ! inout
        &            t_snow_now,     &  ! in
        &            istarts, iends, &  ! in
        &            soiltyp_subs,   &  ! in
        &            w_snow_now      )  ! out

      DO i = istarts, iends
        IF (w_snow_now(i) <= zepsi) THEN
          ! spurious snow is removed
          t_snow_now(i)= t_so_now(i,0)
          w_snow_now(i)= 0.0_ireals
        ELSE
          t_snow_now(i) = MIN( t0_melt, t_snow_now(i))
        ENDIF
      ENDDO

    ENDIF ! multi snow

  ENDIF  ! init_mode

! End of timestep 0 preparations
! ==============================

  END SUBROUTINE terra_multlay_init

  !>
  !! Transform snow depth to snow water equivalent
  !!
  !! Description:
  !!  Snow depth observations are reported in m and this is also the analysis
  !!  variable. The prognostic model variables describing the evolution of
  !!  snow are snow water eqivalent and snow density.
  !!  This subroutine calculates snow water equivalent from snow depth.
  !!  For non-permanent ice points, hsnow is limited to 4m. w_snow is limited
  !!  accordingly, using rho_snow from the first guess.
  !!
  !! Method:
  !!  Simple linear relation between snow water and snow depth using snow density
  !!
  !! @par Revision History
  !! Adapted to ICON by Daniel Reinert, DWD (2014-04-01)
  !! Original implementation by Michael Buchhold (see below)
  !!
  SUBROUTINE get_wsnow(h_snow, rho_snow, t_snow, istart, iend, soiltyp, w_snow)

    REAL(ireals), INTENT(INOUT) :: h_snow(:)    ! snow depth            [m]
    REAL(ireals), INTENT(INOUT) :: rho_snow(:)  ! snow density          [kg/m**3]
    REAL(ireals), INTENT(IN)    :: t_snow(:)    ! snow temperature      [K]
    REAL(ireals), INTENT(INOUT) :: w_snow(:)    ! snow water equivalent [m H2O]
    INTEGER,  INTENT(IN)    :: istart       ! start index
    INTEGER,  INTENT(IN)    :: iend         ! end index
    INTEGER,  INTENT(IN)    :: soiltyp(:)   ! type of soil (1-9)

    ! local
    INTEGER :: jc         ! loop index

    REAL(ireals), PARAMETER :: rho_snw_default=250._ireals
  !--------------------------------------------------------------------

    ! calculate water equivalent from snow depth
    !
    DO jc = istart, iend
      IF (rho_snow(jc) < 1._ireals) rho_snow(jc)=rho_snw_default      ! average initial density
      IF (h_snow(jc) > 0._ireals) THEN
        ! multiply analysed snow depth [m] by (first-guess) density
        ! to get water equivalent in [m H2O]
        IF (soiltyp(jc) /= 1) THEN       ! 1=ice
          ! limit snow depth to 4m for non-glacier points
          h_snow(jc) = MIN(h_snow(jc), 4._ireals)
          w_snow(jc) = h_snow(jc) * rho_snow(jc)/rho_w
        ELSE
          ! limit snow depth to at least 1 m and set snow density to the equilibrium value for the current snow temperature
          h_snow(jc)   = MAX(h_snow(jc), 1._ireals)
          rho_snow(jc) = crhosmax_tmin + MAX(0._ireals,(t_snow(jc)-csnow_tmin)/(t0_melt-csnow_tmin)) &
                         * (crhosmax_ml-crhosmax_tmin)
          w_snow(jc)   = h_snow(jc) * rho_snow(jc)/rho_w
        ENDIF  ! soiltyp
      ENDIF
    ENDDO

  END SUBROUTINE get_wsnow

!!$!
!!$!+ Transform snow depth to snow water equivalent and vice versa
!!$!
!!$      SUBROUTINE snw_tran (snw,rho,ie,je,imode)
!!$!
!!$! Description:
!!$!  Snow depth observations are reported in m and this is also the analysis
!!$!  variable. The prognostic model variables describing the evolution of
!!$!  snow are  snow water eqivalent and snow density.
!!$!  This subroutine calculates snow depth from water eqivalent and vice versa.
!!$!
!!$! Method:
!!$!  Simple linear relation between snow water and snow depth using snow density
!!$!
!!$! Input  files: none
!!$! Output files: none
!!$!
!!$! Current Code Owner: DWD, Michael Buchhold
!!$!  phone:  +49  69  8062 2726
!!$!  fax:    +49  69  8062 3721
!!$!  email:  michael.buchhold@dwd.de
!!$!
!!$! History:
!!$! Version    Date       Name
!!$! ---------- ---------- ----
!!$! 1.4        2003/07/08 Michael Buchhold
!!$!  Initial release
!!$! 1.13       2005/04/18 Michael Buchhold
!!$!  Replacement of common blocks by data modules
!!$! 1.17       2005/12/01 Michael Buchhold
!!$!  Snow density is used for transformation between snow water
!!$!  and snow depth
!!$! V1_26        2010/04/27 Michael Gertz
!!$!  Adaptions to SVN
!!$!
!!$! Code Description:
!!$! Language: Fortran 90.
!!$! Software Standards: "European Standards for Writing and
!!$! Documenting Exchangeable Fortran 90 Code".
!!$!=======================================================================
!!$!
!!$! Declarations:
!!$!
!!$! Modules used:
!!$      USE data_all, ONLY:
!!$     & rhosnw
!!$
!!$      IMPLICIT NONE
!!$!=======================================================================
!!$!
!!$! Subroutine / Function arguments
!!$! Scalar arguments with intent(in):
!!$      INTEGER, INTENT(IN) ::  ie, je   ! Dimensions of snws
!!$      INTEGER, INTENT(IN) ::  imode
!!$!     imode= 1: calculate water equivalent from snow depth
!!$!     imode=-1: calculate snow depth from water equivalent
!!$
!!$! Array arguments with intent(inout):
!!$      REAL, INTENT(INOUT) ::   snw (ie, je), rho(ie,je)
!!$!
!!$! Local Type Definitions:
!!$
!!$!
!!$! Local Parameters:
!!$
!!$!
!!$! Local Scalars:
!!$      INTEGER i,j
!!$      REAL zcrit, zwcrit
!!$!
!!$! Local Arrays:
!!$
!!$!=======================================================================
!!$! Include statements: none
!!$!
!!$!- End of header
!!$!------------------------------------------------------------------------
!!$      if (imode == 1) then
!!$!     calculate water equivalent from snow depth
!!$        do j = 1,je
!!$        do i = 1,ie
!!$          if (snw(i,j) > 0.) then
!!$             if (rho(i,j) < 1.) rho(i,j)=rhosnw
!!$!           scale analysed snow depth [mm] to get in m and multiply it by density
!!$            snw(i,j) = snw(i,j)/1000.*rho(i,j)
!!$          endif
!!$        enddo
!!$        enddo
!!$
!!$      else if (imode ==-1) then
!!$!     calculate  snow depth from water eqivalent
!!$        do j = 1,je
!!$        do i = 1,ie
!!$          if (snw(i,j) >= 0.) then
!!$            if (rho(i,j) < 1.) rho(i,j)=rhosnw
!!$            snw(i,j) = snw(i,j)/rho(i,j)
!!$          else
!!$            rho(i,j) = 0.
!!$          endif
!!$!         scale snow depth to have it in mm
!!$          snw(i,j) = snw(i,j)*1000.
!!$        enddo
!!$        enddo
!!$      endif
!!$
!!$      end

END MODULE mo_nwp_soil_init

max_stars_repo_name: niklaswr/TSMP

max_issues_repo_name: niklaswr/TSMP

max_forks_repo_name: niklaswr/TSMP


text:
Ethnic Restaurants in Davis get a lot of flack for being or not being authentic. This page is intended for discussion of general inauthenticies rather than those specific to restaurants, and can be used as a guide for the style of cuisines a Davis newbie should expect upon arriving here. Restaurant specific criticisms are best suited to individual pages.

Some primary reasons why restaurants seem to be inauthentic:

 1. The original ethnic food might be viewed as undesirable by most Americans and therefore scare away customers. Examples of this include the level of spicyness, use of ingredients not common in American cuisine (e.g.,  Durian fruit, Lutefisk, or Gefilte Fish), organ meats, or meats that resemble the actual animal.
 2. The restaurant skimps on the quality of ingredients, assuming that patrons wont know the difference anyway. Example of this include substituting canned ingredients for fresh ones and adding sugar, MSG, or food coloring to replace expensive flavors rather than to enhance them.
 3. The restaurant might try to cut corners on the proper methods of cooking, to save time and/or money.  The methods used to prepare both ingredients and the dish itself is often very critical in the authenticity of a dish. Perhaps a broth or soup is meant to be simmered for hours, or certain ingredients need to be slowly cooked to bring the right flavor out.  Perhaps they use a food processor rather then finely cutting or mincing it to bring out the right textures in the food.  Maybe the meal is traditionally meant to be cooked over an open flame rather then on stovetop.
 4. The patrons are either unaware of regional variations in the origin area that produced the food or they assume a particular ethnic type when another is being prepared (for instance Southwest versus Texmex versus Mexican... which are each split into dozens of very different subtypes of cuisine).

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!
! NLTECoeff_OLD2NEW
!
! Program to convert netCDF format NLTECoeff files from an OLD release
! format to a NEW one.
!
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 18-Feb-2011
!                       paul.vandelst@noaa.gov
!

PROGRAM NLTECoeff_OLD2NEW

  ! ------------------
  ! Environment set up
  ! ------------------
  ! Module usage
  USE Type_Kinds         , ONLY: fp
  USE File_Utility       , ONLY: Get_Lun, File_Exists
  USE Message_Handler    , ONLY: SUCCESS, FAILURE, WARNING, INFORMATION, &
                                 Display_Message, Program_Message
  USE NLTECoeff_Define   , ONLY: NLTECoeff_type, &
                                 NLTECoeff_Associated, &
                                 NLTECoeff_Create    , &
                                 NLTECoeff_Destroy
  USE NLTECoeff_netCDF_IO, ONLY: NLTECoeff_netCDF_WriteFile
  USE NLTE_Define        , ONLY: NLTE_type, &
                                 NLTE_Create, &
                                 NLTE_Destroy, &
                                 Pt
  USE NLTE_Binary_IO     , ONLY: Read_NLTE_Binary
  USE SensorInfo_Define
  USE SensorInfo_IO
  USE SensorInfo_LinkedList
  ! Disable all implicit typing
  IMPLICIT NONE

  ! ----------
  ! Parameters
  ! ----------
  CHARACTER(*), PARAMETER :: PROGRAM_NAME = 'NLTECoeff_OLD2NEW'
  CHARACTER(*), PARAMETER :: PROGRAM_VERSION_ID = &
  CHARACTER(*), PARAMETER :: SENSORINFO_FILE = 'SensorInfo'

  ! ---------
  ! Variables
  ! ---------
  CHARACTER(256)  :: msg
  CHARACTER(256)  :: old_filename, filename
  INTEGER :: err_stat
  INTEGER :: n
  TYPE(SensorInfo_List_type) :: sinfo_list
  TYPE(SensorInfo_type)      :: sinfo
  TYPE(NLTECoeff_type) :: NLTECoeff
  TYPE(NLTE_type)      :: old_NLTECoeff

  ! Program header
  CALL Program_Message(PROGRAM_NAME, &
                       'Program to convert NLTECoeff files from '//&
                       'an OLD release Binary format to a NEW netCDF one.', &
                       '$Revision$' )

  ! Read the SensorInfo file
  err_stat = Read_SensorInfo( SENSORINFO_FILE, sinfo_list, Quiet = 1 )
  IF ( err_stat /= SUCCESS ) THEN
    msg = 'Error reading SensorInfo file'
    CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
  END IF

  ! Get user inputs
  WRITE( *,FMT='(/5x,"Enter the OLD format NLTECoeff Binary filename: ")',ADVANCE='NO' )
  READ( *,FMT='(a)' ) old_filename
  old_filename = ADJUSTL(old_filename)
  IF ( .NOT. File_Exists( old_filename ) ) THEN
    msg = 'oSRF file '//TRIM(old_filename)//' not found.'
    CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
  END IF
  filename = TRIM(old_filename)//'.nc'

  ! Read the old release netCDF NLTECoeff file
  err_stat = Read_NLTE_Binary( &
    old_filename, &
    old_NLTECoeff )
  IF ( err_stat /= SUCCESS ) THEN
    msg = 'Error reading old release binary NLTECoeff file '//TRIM(old_filename)
    CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
  END IF

  ! Get the matching SensorInfo list node
  Node_Loop: DO n = 1, Count_SensorInfo_Nodes( sinfo_list )
    err_stat = GetFrom_SensorInfo_List( sinfo_list, n, sinfo )
    IF ( err_stat /= SUCCESS ) THEN
      WRITE( msg,'("Error retrieving node #",i0," from SensorInfo file")' ) n
      CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
    END IF
    IF ( TRIM(sinfo%Sensor_Id) == TRIM(old_NLTECoeff%Sensor_Id) ) EXIT Node_Loop
    err_stat = Destroy_SensorInfo( sinfo )
  END DO Node_Loop
  ! ...Check node is valid
  IF ( .NOT. Associated_SensorInfo( sinfo ) ) THEN
    msg = 'Error finding SensorInfo node for '//TRIM(old_NLTECoeff%Sensor_Id)
    CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
  END IF

  ! Allocate the new NLTECoeff structure
  CALL NLTECoeff_Create( &
    NLTECoeff, &
    old_NLTECoeff%n_Predictors    , &
    old_NLTECoeff%n_Sensor_ZAngles, &
    old_NLTECoeff%n_Solar_ZAngles , &
    old_NLTECoeff%n_NLTE_Channels , &
    old_NLTECoeff%n_Channels        )
  IF ( .NOT. NLTECoeff_Associated( NLTECoeff ) ) THEN
    msg = 'Error allocating new NLTECoeff strructure'
    CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
  END IF

  ! Copy over the unchanged data
  NLTECoeff%Sensor_Id            = sinfo%Sensor_Id
  NLTECoeff%WMO_Satellite_ID     = sinfo%WMO_Satellite_ID
  NLTECoeff%WMO_Sensor_ID        = sinfo%WMO_Sensor_ID
  NLTECoeff%Sensor_Channel       = sinfo%Sensor_Channel
  NLTECoeff%Upper_Plevel         = Pt(1,:)
  NLTECoeff%Lower_Plevel         = Pt(2,:)
  NLTECoeff%Min_Tm               = old_NLTECoeff%Tm(1,:)
  NLTECoeff%Max_Tm               = old_NLTECoeff%Tm(2,:)
  NLTECoeff%Mean_Tm              = old_NLTECoeff%Tm(3,:)
  NLTECoeff%C_Index              = old_NLTECoeff%NLTE_Channel_Index
  NLTECoeff%C                    = old_NLTECoeff%C
  NLTECoeff%Secant_Sensor_Zenith = old_NLTECoeff%Sensor_Zangle
  NLTECoeff%Secant_Solar_Zenith  = old_NLTECoeff%Solar_Zangle
  NLTECoeff%NLTE_Channel         = old_NLTECoeff%NLTE_Channel
  WHERE(old_NLTECoeff%NLTE_Channel_Index > 0) NLTECoeff%Is_NLTE_Channel = .TRUE.

  ! Write the new netCDF NLTECoeff file
  err_stat = NLTECoeff_netCDF_WriteFile( &
    filename, &
    NLTECoeff, &
    Title   = 'NLTE-correction coefficients for '//TRIM(NLTECoeff%Sensor_Id), &
    History = PROGRAM_VERSION_ID, &
    Comment = 'Created from old version datafiles for testing.' )
  IF ( err_stat /= SUCCESS ) THEN
    msg = 'Error writing new release netCDF NLTECoeff file '//TRIM(filename)
    CALL Display_Message( PROGRAM_NAME, msg, FAILURE ); STOP
  END IF

  ! Clean up
  err_stat = Destroy_SensorInfo( sinfo )
  err_stat = Destroy_SensorInfo_List( sinfo_list, Quiet = 1 )
  CALL NLTECoeff_Destroy( NLTECoeff )
  CALL NLTE_Destroy( old_NLTECoeff )

END PROGRAM NLTECoeff_OLD2NEW

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


text:
C****
C                    ************************
C                    *     OASIS MODULE     *
C                    *     ------------     *
C                    ************************
C****
C***********************************************************************
C     This module belongs to the SCRIP library. It is a modified version
C     of SCRIP 1.4 remap_ditswgt.f in order to weight a neighbour by
C     exp[-1/2 * d^^2/sigma^^2] where d is the distance between the source
C     and the target points and sigma^^2=VAR*dm^^2, where dm is the
C     average distance between the source grid points and VAR is a value
C     given by the user.
C
C     Additional Modifications:
C       - restrict types are written in capital letters
C       - bug line 429: bin_lons(3,n) instead of bin_lons(2,n)
C
C     Modified by            D. Declat,  CERFACS              27.06.2002
C***********************************************************************
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!     this module contains necessary routines for performing an
!     interpolation using a distance-weighted average of n nearest
!     neighbors.
!
!-----------------------------------------------------------------------
!
!     CVS:$Id: remap_gauswgt.f 2826 2010-12-10 11:14:21Z valcke $
!
!     Copyright (c) 1997, 1998 the Regents of the University of
!       California.
!
!     This software and ancillary information (herein called software)
!     called SCRIP is made available under the terms described here.
!     The software has been approved for release with associated
!     LA-CC Number 98-45.
!
!     Unless otherwise indicated, this software has been authored
!     by an employee or employees of the University of California,
!     operator of the Los Alamos National Laboratory under Contract
!     No. W-7405-ENG-36 with the U.S. Department of Energy.  The U.S.
!     Government has rights to use, reproduce, and distribute this
!     software.  The public may copy and use this software without
!     charge, provided that this Notice and any statement of authorship
!     are reproduced on all copies.  Neither the Government nor the
!     University makes any warranty, express or implied, or assumes
!     any liability or responsibility for the use of this software.
!
!     If software is modified to produce derivative works, such modified
!     software should be clearly marked, so as not to confuse it with
!     the version available from Los Alamos National Laboratory.
!
!***********************************************************************

      module remap_gaussian_weight

!-----------------------------------------------------------------------

      use kinds_mod     ! defines common data types
      use constants     ! defines common constants
      use grids         ! module containing grid info
      use remap_vars    ! module containing remap info
      USE mod_oasis_flush

      implicit none

!-----------------------------------------------------------------------
!
!     module variables
!
!-----------------------------------------------------------------------

      real (kind=dbl_kind), dimension(:), allocatable, save ::
     &     coslat, sinlat, ! cosine, sine of grid lats (for distance)
     &     coslon, sinlon, ! cosine, sine of grid lons (for distance)
     &     wgtstmp         ! an array to hold the link weight

!***********************************************************************

      contains

!***********************************************************************

      subroutine remap_gauswgt (lextrapdone, num_neighbors, rl_varmul)

!-----------------------------------------------------------------------
!
!     this routine computes the inverse-distance weights for a
!     nearest-neighbor interpolation.
!
!-----------------------------------------------------------------------

      LOGICAL ::
     &           lextrapdone   ! logical, true if EXTRAP done on field

      REAL (kind=real_kind) ::
     $    rl_varmul             ! Gaussian variance

      INTEGER (kind=int_kind) ::
     &       num_neighbors     ! number of neighbours

!-----------------------------------------------------------------------
!
!     local variables
!
!-----------------------------------------------------------------------

      logical (kind=log_kind), dimension(num_neighbors) ::
     &     nbr_mask        ! mask at nearest neighbors

      integer (kind=int_kind) :: n,
     &     dst_add,        ! destination address
     &     nmap            ! index of current map being computed

      integer (kind=int_kind), dimension(num_neighbors) ::
     &     nbr_add         ! source address at nearest neighbors

      real (kind=dbl_kind), dimension(num_neighbors) ::
     &     nbr_dist        ! angular distance four nearest neighbors

      real (kind=dbl_kind) ::
     &     coslat_dst,     ! cos(lat) of destination grid point
     &     coslon_dst,     ! cos(lon) of destination grid point
     &     sinlat_dst,     ! sin(lat) of destination grid point
     &     sinlon_dst,     ! sin(lon) of destination grid point
     &     dist_tot,       ! sum of neighbor distances (for normalizing)
     &     dist_average    ! average distance between the source points
      logical (kind=log_kind) :: ll_allmask
      logical (kind=log_kind), PARAMETER :: ll_nnei=.true.
      real (kind=dbl_kind) ::
     &        distance ,plat,plon,src_latsnn, arg   ! angular distance
      real (kind=dbl_kind), dimension (1) :: wgts_new
      integer (kind=int_kind) :: min_add_out, max_add_out,
     &        srch_add, src_addnn
!-----------------------------------------------------------------------
!
      IF (nlogprt .GE. 2) THEN
         WRITE (UNIT = nulou,FMT = *)' '
         WRITE (UNIT = nulou,FMT = *)'Entering routine remap_gauswgt'
         CALL OASIS_FLUSH_SCRIP(nulou)
      ENDIF
!-----------------------------------------------------------------------
!
!     compute mappings from grid1 to grid2
!
!-----------------------------------------------------------------------

      nmap = 1

      !***
      !*** allocate wgtstmp to be consistent with store_link interface
      !***

      allocate (wgtstmp(num_wts))

      !***
      !*** compute cos, sin of lat/lon on source grid for distance
      !*** calculations
      !***

      allocate (coslat(grid1_size), coslon(grid1_size),
     &          sinlat(grid1_size), sinlon(grid1_size))

      coslat = cos(grid1_center_lat)
      coslon = cos(grid1_center_lon)
      sinlat = sin(grid1_center_lat)
      sinlon = sin(grid1_center_lon)

      !***
      !*** compute the average of the distances between the source points
      !***

      call grid_dist_average(grid1_size,
     &                       coslat, coslon,
     &                       sinlat, sinlon,
     &                       dist_average)
      !***
      !*** loop over destination grid
      !***

      grid_loop1: do dst_add = 1, grid2_size

        if (.not. grid2_mask(dst_add)) cycle grid_loop1

        coslat_dst = cos(grid2_center_lat(dst_add))
        coslon_dst = cos(grid2_center_lon(dst_add))
        sinlat_dst = sin(grid2_center_lat(dst_add))
        sinlon_dst = sin(grid2_center_lon(dst_add))

        !***
        !*** find nearest grid points on source grid and
        !*** distances to each point
        !***

        call grid_search_nbrg(nbr_add, nbr_dist,
     &                        min_add_out, max_add_out,
     &                        grid2_center_lat(dst_add),
     &                        grid2_center_lon(dst_add),
     &                        coslat_dst, coslon_dst,
     &                        sinlat_dst, sinlon_dst,
     &                        bin_addr1,
     &                        dist_average,  num_neighbors, rl_varmul)

        !***
        !*** compute weights based on inverse distance
        !*** if mask is false, eliminate those points
        !***
        src_addnn = zero
        dist_tot = zero
        do n=1,num_neighbors
          if ((grid1_mask(nbr_add(n))) .or.
     &        (.not. grid1_mask(nbr_add(n)) .and. lextrapdone)) THEN
            nbr_dist(n) = one/nbr_dist(n)
            dist_tot = dist_tot + nbr_dist(n)
            nbr_mask(n) = .true.
          else
            nbr_mask(n) = .false.
          endif
        end do

        !***
        !*** normalize weights and store the link
        !***

        do n=1,num_neighbors
          if (nbr_mask(n)) then
            wgtstmp(1) = nbr_dist(n)/dist_tot
            call store_link_nbr(nbr_add(n), dst_add, wgtstmp, nmap)
            grid2_frac(dst_add) = one
          endif
        end do
        IF (ll_nnei) THEN
            ll_allmask= .true.
            do n=1,num_neighbors
              if (nbr_mask(n)) then
                  ll_allmask=.false.
              endif
            end do
            if ( ll_allmask) THEN
                src_latsnn = bignum
                do srch_add = min_add_out,max_add_out
                  if (grid1_mask(srch_add) .or.
     &                (.not. grid1_mask(srch_add)
     &                .and. lextrapdone)) THEN
                      arg = coslat_dst*cos(grid1_center_lat(srch_add))*
     &                (coslon_dst*cos(grid1_center_lon(srch_add)) +
     &                sinlon_dst*sin(grid1_center_lon(srch_add)))+
     &                sinlat_dst*sin(grid1_center_lat(srch_add))
                      IF (arg < -1.0d0) THEN
                          arg = -1.0d0
                      ELSE IF (arg > 1.0d0) THEN
                          arg = 1.0d0
                      END IF
                      distance=acos(arg)

                      if (distance < src_latsnn) then
                          src_addnn = srch_add
                          src_latsnn = distance
                      endif
                  endif
                end do
                IF (nlogprt .GE. 2) THEN
                    WRITE(nulou,*)'ll_allmask=true and src_addnn= '
     &                  ,src_addnn
                    CALL OASIS_FLUSH_SCRIP(nulou)
                ENDIF
                wgts_new(1) = 1.
                grid2_frac(dst_add) = one
                call store_link_nbr(src_addnn,dst_add ,wgts_new, nmap)
            endif
        endif

      end do grid_loop1

      deallocate (coslat, coslon, sinlat, sinlon)

!-----------------------------------------------------------------------
!
!     compute mappings from grid2 to grid1 if necessary
!
!-----------------------------------------------------------------------

      if (num_maps > 1) then

      nmap = 2

      !***
      !*** compute cos, sin of lat/lon on source grid for distance
      !*** calculations
      !***

      allocate (coslat(grid2_size), coslon(grid2_size),
     &          sinlat(grid2_size), sinlon(grid2_size))

      coslat = cos(grid2_center_lat)
      coslon = cos(grid2_center_lon)
      sinlat = sin(grid2_center_lat)
      sinlon = sin(grid2_center_lon)

      !***
      !*** compute the average of the distances between the source points
      !***

      call grid_dist_average(grid2_size,
     &                       coslat, coslon,
     &                       sinlat, sinlon,
     &                       dist_average)

      !***
      !*** loop over destination grid
      !***

      grid_loop2: do dst_add = 1, grid1_size

        if (.not. grid1_mask(dst_add)) cycle grid_loop2

        coslat_dst = cos(grid1_center_lat(dst_add))
        coslon_dst = cos(grid1_center_lon(dst_add))
        sinlat_dst = sin(grid1_center_lat(dst_add))
        sinlon_dst = sin(grid1_center_lon(dst_add))

        !***
        !*** find four nearest grid points on source grid and
        !*** distances to each point
        !***

        call grid_search_nbrg(nbr_add, nbr_dist,
     &                        min_add_out, max_add_out,
     &                        grid1_center_lat(dst_add),
     &                        grid1_center_lon(dst_add),
     &                        coslat_dst, coslon_dst,
     &                        sinlat_dst, sinlon_dst,
     &                        bin_addr2,
     &                        dist_average, num_neighbors, rl_varmul)

        !***
        !*** compute weights based on inverse distance
        !*** if mask is false, eliminate those points
        !***

        dist_tot = zero
        do n=1,num_neighbors
          if (grid2_mask(nbr_add(n))) then
            nbr_dist(n) = one/nbr_dist(n)
            dist_tot = dist_tot + nbr_dist(n)
            nbr_mask(n) = .true.
          else
            nbr_mask(n) = .false.
          endif
        end do

        !***
        !*** normalize weights and store the link
        !***

        do n=1,num_neighbors
          if (nbr_mask(n)) then
            wgtstmp(1) = nbr_dist(n)/dist_tot
            call store_link_nbr(dst_add, nbr_add(n), wgtstmp, nmap)
            grid1_frac(dst_add) = one
          endif
        end do

        IF (ll_nnei) THEN
            ll_allmask= .true.
            do n=1,num_neighbors
              if (nbr_mask(n)) then
                  ll_allmask=.false.
              endif
            end do
            if ( ll_allmask) then
                PRINT*, 'll_allmask true',src_addnn
                src_latsnn = bignum
                do srch_add = min_add_out,max_add_out
                  if (grid2_mask(srch_add) .or.
     &                (.not. grid2_mask(srch_add)
     &                .and. lextrapdone)) THEN
                      arg = coslat_dst*cos(grid2_center_lat(srch_add))*
     &                (coslon_dst*cos(grid2_center_lon(srch_add)) +
     &                sinlon_dst*sin(grid2_center_lon(srch_add)))+
     &                sinlat_dst*sin(grid2_center_lat(srch_add))
                      IF (arg < -1.0d0) THEN
                          arg = -1.0d0
                      ELSE IF (arg > 1.0d0) THEN
                          arg = 1.0d0
                      END IF
                      distance=acos(arg)

                      if (distance < src_latsnn) then
                          src_addnn = srch_add
                          src_latsnn = distance
                      endif
                  endif
                end do
                wgts_new = 1.
                grid1_frac(dst_add) = one
                call store_link_nbr(dst_add, src_addnn, wgts_new, nmap)
            endif
        endif
      end do grid_loop2

      deallocate (coslat, coslon, sinlat, sinlon)

      endif

      deallocate(wgtstmp)
!
      IF (nlogprt .GE. 2) THEN
         WRITE (UNIT = nulou,FMT = *)' '
         WRITE (UNIT = nulou,FMT = *)'Leaving routine remap_gauswgt'
         CALL OASIS_FLUSH_SCRIP(nulou)
      ENDIF
!
!-----------------------------------------------------------------------

      end subroutine remap_gauswgt

!***********************************************************************

      subroutine grid_search_nbrg(nbr_add, nbr_dist, min_add, max_add,
     &               plat, plon,
     &               coslat_dst, coslon_dst, sinlat_dst, sinlon_dst,
     &               src_bin_add, dst_average,
     $               num_neighbors, rl_varmul)

!-----------------------------------------------------------------------
!
!     this routine finds the closest num_neighbor points to a search
!     point and computes a distance to each of the neighbors.
!
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!
!     input variables
!
!-----------------------------------------------------------------------

      integer (kind=int_kind), dimension(:,:), intent(in) ::
     &        src_bin_add ! search bins for restricting search

      real (kind=dbl_kind), intent(in) ::
     &        plat,         ! latitude  of the search point
     &        plon,         ! longitude of the search point
     &        coslat_dst,   ! cos(lat)  of the search point
     &        coslon_dst,   ! cos(lon)  of the search point
     &        sinlat_dst,   ! sin(lat)  of the search point
     &        sinlon_dst,   ! sin(lon)  of the search point
     &        dst_average   ! average distance between the source points

      REAL (kind=real_kind), intent(in) ::
     &        rl_varmul     ! Gaussian variance

      INTEGER (kind=int_kind) ::
     &       num_neighbors     ! number of neighbours

!-----------------------------------------------------------------------
!
!     output variables
!
!-----------------------------------------------------------------------

      integer (kind=int_kind), dimension(num_neighbors), intent(out) ::
     &        nbr_add  ! address of each of the closest points

      real (kind=dbl_kind), dimension(num_neighbors), intent(out) ::
     &        nbr_dist ! distance to each of the closest points

      integer (kind=int_kind),intent(out) :: min_add, max_add

!-----------------------------------------------------------------------
!
!     local variables
!
!-----------------------------------------------------------------------

      integer (kind=int_kind) :: n, nmax, nadd, nchk, ! dummy indices
     &        nm1, np1, i, j, ip1, im1, jp1, jm1

      real (kind=dbl_kind) ::
     &        distance, arg      ! angular distance

      real (kind=dbl_kind) ::
     &        variance      ! variance for the gaussian FUNCTION

!-----------------------------------------------------------------------
!
!     loop over source grid and find nearest neighbors
!
!-----------------------------------------------------------------------

      !***
      !*** restrict the search using search bins
      !*** expand the bins to catch neighbors
      !***

      select case (restrict_type)
      case('LATITUDE')

        do n=1,num_srch_bins
          if (plat >= bin_lats(1,n) .and. plat <= bin_lats(2,n)) then
            min_add = src_bin_add(1,n)
            max_add = src_bin_add(2,n)

            nm1 = max(n-1,1)
            np1 = min(n+1,num_srch_bins)

            min_add = min(min_add,src_bin_add(1,nm1))
            max_add = max(max_add,src_bin_add(2,nm1))
            min_add = min(min_add,src_bin_add(1,np1))
            max_add = max(max_add,src_bin_add(2,np1))
          endif
        end do

      case('LATLON')

        n = 0
        nmax = nint(sqrt(real(num_srch_bins)))
        do j=1,nmax
        jp1 = min(j+1,nmax)
        jm1 = max(j-1,1)
        do i=1,nmax
          ip1 = min(i+1,nmax)
          im1 = max(i-1,1)

          n = n+1
          if (plat >= bin_lats(1,n) .and. plat <= bin_lats(2,n) .and.
     &        plon >= bin_lons(1,n) .and. plon <= bin_lons(2,n)) then
            min_add = src_bin_add(1,n)
            max_add = src_bin_add(2,n)

            nm1 = (jm1-1)*nmax + im1
            np1 = (jp1-1)*nmax + ip1
            nm1 = max(nm1,1)
            np1 = min(np1,num_srch_bins)

            min_add = min(min_add,src_bin_add(1,nm1))
            max_add = max(max_add,src_bin_add(2,nm1))
            min_add = min(min_add,src_bin_add(1,np1))
            max_add = max(max_add,src_bin_add(2,np1))
          endif
        end do
        end do

      end select

      !***
      !*** initialize distance and address arrays
      !***

      nbr_add = 0
      nbr_dist = bignum
      variance = rl_varmul*dst_average*dst_average
      do nadd=min_add,max_add

        !***
        !*** find distance to this point
        !***
        arg = sinlat_dst*sinlat(nadd) +
     &        coslat_dst*coslat(nadd)*
     &       (coslon_dst*coslon(nadd) +
     &        sinlon_dst*sinlon(nadd))
        IF (arg < -1.0d0) THEN
            arg = -1.0d0
        ELSE IF (arg > 1.0d0) THEN
            arg = 1.0d0
        END IF
        distance = acos(arg)
        !distance = min(500., distance) !ts-sv
        !distance = exp(.5*distance*distance/variance)
        !***
        !*** store the address and distance if this is one of the
        !*** smallest four so far
        !***

        check_loop: do nchk=1,num_neighbors
          if (distance .lt. nbr_dist(nchk)) THEN
            do n=num_neighbors,nchk+1,-1
              nbr_add(n) = nbr_add(n-1)
              nbr_dist(n) = nbr_dist(n-1)
            end do
            nbr_add(nchk) = nadd
            nbr_dist(nchk) = distance
            exit check_loop
          endif
        end do check_loop

      end do

      exp_loop: do nchk=1,num_neighbors
          nbr_dist(nchk) =
     &     exp(.5*nbr_dist(nchk)*nbr_dist(nchk)/variance)
      end do exp_loop

!-----------------------------------------------------------------------

      end subroutine grid_search_nbrg

!***********************************************************************

      subroutine store_link_nbr(add1, add2, weights, nmap)

!-----------------------------------------------------------------------
!
!     this routine stores the address and weight for this link in
!     the appropriate address and weight arrays and resizes those
!     arrays if necessary.
!
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!
!     input variables
!
!-----------------------------------------------------------------------

      integer (kind=int_kind), intent(in) ::
     &        add1,  ! address on grid1
     &        add2,  ! address on grid2
     &        nmap   ! identifies which direction for mapping

      real (kind=dbl_kind), dimension(:), intent(in) ::
     &        weights ! array of remapping weights for this link

!-----------------------------------------------------------------------
!
!     increment number of links and check to see if remap arrays need
!     to be increased to accomodate the new link.  then store the
!     link.
!
!-----------------------------------------------------------------------

      select case (nmap)
      case(1)

        num_links_map1  = num_links_map1 + 1

        if (num_links_map1 > max_links_map1)
     &     call resize_remap_vars(1,resize_increment)

        grid1_add_map1(num_links_map1) = add1
        grid2_add_map1(num_links_map1) = add2
        wts_map1    (:,num_links_map1) = weights

      case(2)

        num_links_map2  = num_links_map2 + 1

        if (num_links_map2 > max_links_map2)
     &     call resize_remap_vars(2,resize_increment)

        grid1_add_map2(num_links_map2) = add1
        grid2_add_map2(num_links_map2) = add2
        wts_map2    (:,num_links_map2) = weights

      end select

!-----------------------------------------------------------------------

      end subroutine store_link_nbr

!***********************************************************************

      subroutine grid_dist_average(grid_size,
     &                             coslat_grid, coslon_grid,
     &                             sinlat_grid, sinlon_grid,
     &                             dist_average)

!-----------------------------------------------------------------------
!
!     this routine computes the average distance between the points of a
!     grid.
!
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!
!     output variables
!
!-----------------------------------------------------------------------

      real (kind=dbl_kind), intent(out) ::
     &        dist_average ! distance to each of the closest points

!-----------------------------------------------------------------------
!
!     input variables
!
!-----------------------------------------------------------------------

      integer (kind=int_kind), intent(in) ::
     &        grid_size

      real (kind=dbl_kind), dimension(:), intent(in) ::
     &        coslat_grid,   ! cos(lat)  of the grid points
     &        coslon_grid,   ! cos(lon)  of the grid points
     &        sinlat_grid,   ! sin(lat)  of the grid points
     &        sinlon_grid    ! sin(lon)  of the grid points
      REAL (kind=dbl_kind) :: arg

!-----------------------------------------------------------------------
!
!     local variables
!
!-----------------------------------------------------------------------

      integer (kind=int_kind) :: i
!
      IF (nlogprt .GE. 2) THEN
         WRITE (UNIT = nulou,FMT = *)' '
         WRITE (UNIT = nulou,FMT = *)
     &       'Entering routine remap_dist_average'
         CALL OASIS_FLUSH_SCRIP(nulou)
      ENDIF
!
!-----------------------------------------------------------------------
!
!     compute the distance over the grid and average
!
!-----------------------------------------------------------------------

      dist_average = 0.0
      DO i = 1, grid_size
        IF (i .eq. 1) THEN
            arg = sinlat_grid(grid_size)*sinlat_grid(i) +
     &          coslat_grid(grid_size)*coslat_grid(i)*
     &          (coslon_grid(grid_size)*coslon_grid(i) +
     &          sinlon_grid(grid_size)*sinlon_grid(i))
            IF (arg < -1.0d0) THEN
                arg = -1.0d0
            ELSE IF (arg > 1.0d0) THEN
                arg = 1.0d0
            END IF
            dist_average = dist_average + acos(arg)

            arg = sinlat_grid(i)*sinlat_grid(i+1) +
     &          coslat_grid(i)*coslat_grid(i+1)*
     &          (coslon_grid(i)*coslon_grid(i+1) +
     &          sinlon_grid(i)*sinlon_grid(i+1))
            IF (arg < -1.0d0) THEN
                arg = -1.0d0
            ELSE IF (arg > 1.0d0) THEN
                arg = 1.0d0
            END IF
            dist_average = dist_average + acos(arg)

        ELSE IF (i .eq. grid_size) THEN
            arg = sinlat_grid(i-1)*sinlat_grid(i) +
     &          coslat_grid(i-1)*coslat_grid(i)*
     &          (coslon_grid(i-1)*coslon_grid(i) +
     &          sinlon_grid(i-1)*sinlon_grid(i))
            IF (arg < -1.0d0) THEN
                arg = -1.0d0
            ELSE IF (arg > 1.0d0) THEN
                arg = 1.0d0
            END IF
            dist_average = dist_average + acos(arg)

            arg = sinlat_grid(i)*sinlat_grid(1) +
     &          coslat_grid(i)*coslat_grid(1)*
     &          (coslon_grid(i)*coslon_grid(1) +
     &          sinlon_grid(i)*sinlon_grid(1))
            IF (arg < -1.0d0) THEN
                arg = -1.0d0
            ELSE IF (arg > 1.0d0) THEN
                arg = 1.0d0
            END IF
            dist_average = dist_average + acos(arg)

        ELSE
            arg = sinlat_grid(i-1)*sinlat_grid(i) +
     &          coslat_grid(i-1)*coslat_grid(i)*
     &          (coslon_grid(i-1)*coslon_grid(i) +
     &          sinlon_grid(i-1)*sinlon_grid(i))
            IF (arg < -1.0d0) THEN
                arg = -1.0d0
            ELSE IF (arg > 1.0d0) THEN
                arg = 1.0d0
            END IF
            dist_average = dist_average + acos(arg)

            arg = sinlat_grid(i)*sinlat_grid(i+1) +
     &          coslat_grid(i)*coslat_grid(i+1)*
     &          (coslon_grid(i)*coslon_grid(i+1) +
     &          sinlon_grid(i)*sinlon_grid(i+1))
            IF (arg < -1.0d0) THEN
                arg = -1.0d0
            ELSE IF (arg > 1.0d0) THEN
                arg = 1.0d0
            END IF
            dist_average = dist_average + acos(arg)
        END IF
      END DO
      dist_average = dist_average / (2*grid_size)
!
      IF (nlogprt .GE. 2) THEN
         WRITE (UNIT = nulou,FMT = *)' '
         WRITE (UNIT = nulou,FMT = *)
     &       'Leaving routine remap_dist_average'
         CALL OASIS_FLUSH_SCRIP(nulou)
      ENDIF
!
      END subroutine grid_dist_average

!***********************************************************************

      end module remap_gaussian_weight

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

max_stars_repo_name: minsukji/ci-debug

max_issues_repo_name: minsukji/ci-debug

max_forks_repo_name: minsukji/ci-debug


module pinfo_module

  implicit none
  private
  public :: get_world_rank_pinfo

  integer,allocatable :: node_partition(:)
  integer :: irank,irank_world
  logical :: init_done=.false.
  integer :: icheck(4)

contains

  subroutine init
    use io_tools_module, only: IOTools_parseInteger
    implicit none
    call IOTools_parseInteger( 'PROCS', node_partition )
    init_done = .true.
  end subroutine init

  integer function get_world_rank_pinfo( g, b, k, s )
    implicit none
    integer,intent(in) :: g, b, k, s
    integer,allocatable :: itmp(:)
    integer :: n
    if ( .not.init_done ) call init
    n=size(node_partition)
    allocate( itmp(n) ); itmp=0
    irank=-1
    irank_world=-1
    icheck = (/ g, b, k, s/)
    call do_loop( n, node_partition, itmp )
    get_world_rank_pinfo = irank_world
  end function get_world_rank_pinfo

  recursive subroutine do_loop( n, nnn, itmp )
    implicit none
    integer,intent(in) :: n
    integer,intent(in) :: nnn(:)
    integer,intent(inout) :: itmp(:)
    integer :: i,irank_gbks(4)
    do i = 0, nnn(n)-1
      itmp(n)=i
      if ( n > 1 ) then
        call do_loop( n-1, nnn, itmp )
      else
        irank = irank + 1
        irank_gbks(1) = itmp(1) + itmp(2)*nnn(1) + itmp(3)*product(nnn(1:2))
        irank_gbks(2:4) = itmp(4:6)
        if ( all(irank_gbks == icheck) ) irank_world = irank
      end if
    end do
  end subroutine do_loop

end module pinfo_module


text:
      SUBROUTINE SINQI (N,WSAVE)
      DIMENSION       WSAVE(1)
      CALL COSQI (N,WSAVE)
      RETURN
      END

max_stars_repo_name: mreineck/benchfft

max_issues_repo_name: mreineck/benchfft

max_forks_repo_name: mreineck/benchfft


text:
program main
	use sky_util
	implicit none
	real t0,t1,temp
	real angle
	character arg*20
!	real ref
!	real,dimension(imax,jmax) :: REF

	if (iargc().ne.2) stop 'Argument number should be 3 !'
	call getarg(1,arg)
!	open(1,file="kiban50_450F.futm",form='unformatted')
	open(1,file=arg,form='unformatted')
	read(1) DEM
	close(1)
	call getarg(2,arg)
	read(arg,*) angle
	call cpu_time(t0)
	call view(angle)
	call cpu_time(temp)
	write(*,*) temp-t0,'sec'
	!open(2,file="dem2.img",form='unformatted')
	!write(2) DEM2
	!close(2)
	!open(3,file="view2.img",form='unformatted')
	!write(3) DVIEW2
	!close(3)
	open(4,file="angle.img",form='unformatted')
	write(4) DVIEW
	close(4)

end program main

max_stars_repo_name: y-iikura/TerrainParam

max_issues_repo_name: y-iikura/TerrainParam

max_forks_repo_name: y-iikura/TerrainParam