Subroutine vom_get_soilprofile

683 : 
684 : !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
685 : !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
686 : !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
687 : !*-----PARAMETER READING FROM SOILPROFILE.PAR---------------------------
688 : 
689 :       subroutine vom_get_soilprofile ()
690 :       use vom_vegwat_mod
691 :       implicit none
692 : 
693 :       INTEGER :: iostat, i, j, indlayer
694 :       LOGICAL :: foundlayer
695 : 
696 : !     * The file soilprofile.par can contain information about thickness
697 : !       and soil properties in each soil layer, with the layer number in
698 : !       the first column.
699 : 
700 :       open(kfile_soilprofile, FILE=trim(adjustl(i_inputpath))// &
701 :            trim(adjustl(sfile_soilprofile)),                  &
702 :      &                        STATUS='old', IOSTAT=iostat)
703 :       if (iostat .eq. 0) then
704 :         do j = 1, s_maxlayer
705 :           read(kfile_soilprofile,*) s_maxlayer, s_delz(j), s_ksat(j),  &
706 :      &      s_nvg(j), s_avg(j), s_thetas(j), s_thetar(j)
707 :         enddo
708 : 
709 : !        Check if i_cz aligns with s_delz
710 : !        Raise a warning and correct if this is not the case
711 :          if ( abs( nint(i_cz / i_delz) - (i_cz / i_delz) )  .gt. 1.0d-6) then
712 :            write(*,*) "ERROR: i_cz does not align with soil layers"
713 :            write(*,*) " Please correct in soilprofile.par and restart"
714 :            stop
715 :          end if
716 : 
717 : !        Check if i_zr aligns with s_delz
718 : !        Raise a warning and correct if this is not the case
719 : 
720 : !       Find the layer of i_zr
721 :         foundlayer = .FALSE.
722 :         do j = 1, s_maxlayer
723 :            if( (abs(sum(s_delz(1:j)) - (i_cz - i_zr) ) .lt. 1.0d-6) .and. &
724 :                (foundlayer .eqv. .FALSE.)  ) then
725 :                indlayer = j
726 :                foundlayer = .TRUE.
727 :            end if
728 :         enddo
729 : 
730 : !       correting i_zr
731 :         if( abs(sum(s_delz(1:indlayer)) - (i_cz - i_zr) ) .gt. 1.0d-6 ) then
732 :            write(*,*) "ERROR: i_zr does not align with soil layers"
733 :            write(*,*) " Please correct in soilprofile.par and restart"
734 :          stop
735 :         end if
736 : 
737 :         c_mvg(:) = 1.d0 - (1.d0 / s_nvg(:))  ! van Genuchten soil parameter m
738 :       else
739 :         s_delz(:)   = i_delz
740 :         s_ksat(:)   = i_ksat
741 :         s_nvg(:)    = i_nvg
742 :         s_avg(:)    = i_avg
743 :         s_thetas(:) = i_thetas
744 :         s_thetar(:) = i_thetar
745 :         c_mvg(:)    = i_mvg
746 :       endif
747 :       close(kfile_soilprofile)
748 : 
749 :       do i = 1, s_maxlayer
750 :         c_hhydrst(i) = (i - 0.5d0) * s_delz(i)  ! (Out[238]) hydrostatic head for (3.34)
751 :       enddo
752 : 
753 :       return
754 :       end subroutine vom_get_soilprofile