Subroutine vom_calc_derived

889 : 
890 : !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
891 : !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
892 : !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
893 : !*-----Calculation of derived parameters--------------------------------
894 : 
895 :       subroutine vom_calc_derived ()
896 :       use vom_vegwat_mod
897 :       implicit none
898 : 
899 :       INTEGER :: in, ik, ii
900 :       INTEGER :: in1, in2
901 :       REAL*8  :: sunr, suns
902 :       REAL*8  :: tairmean
903 :       REAL*8  :: dtair
904 :       REAL*8  :: daylength              ! Day length (hours)
905 :       REAL*8  :: vp__                   ! Absolute vapour pressure in the air (Pa)
906 : 
907 :         in1 = 1
908 :         in2 = c_maxday
909 : 
910 :       if (i_write_h == 1) then
911 :         open(kfile_hourlyweather, FILE=trim(adjustl(i_inputpath))// &
912 :              trim(adjustl(sfile_hourlyweather)), STATUS='new')
913 :         write(kfile_hourlyweather,'(5a8,5a11)') 'hour', 'dayyear', 'fday', &
914 :      &    'fmonth', 'fyear', 'tair_h', 'vd_h', 'par_h', 'rain_h', 'ca_h'
915 :       endif
916 : 
917 :       do in = in1, in2
918 :         par_d(in) = 2.0804d0 * srad_d(in)  ! (Out[17]), par in mol/m2 if srad was MJ/m2
919 :         daylength = 12.d0 - 7.639437d0 * ASIN((0.397949d0              &
920 :      &            * COS(0.172142d0 + 0.017214d0 * dayyear(in))         &
921 :      &            * TAN(0.017453d0 * i_lat))                           &
922 :      &            / SQRT(0.920818d0 - 0.079182d0                       &
923 :      &            * COS(0.344284d0 + 0.034428d0 * dayyear(in))))  ! (Out[22]), in hours
924 : !       * sets time of sunrise and sunset
925 :         sunr = 12d0 - 0.5d0 * daylength
926 :         suns = 12d0 + 0.5d0 * daylength
927 :         tairmean = (tairmax_d(in) + tairmin_d(in)) / 2.d0
928 :         dtair = tairmax_d(in) - tairmin_d(in)
929 :         vp__ = vp_d(in) * 100.d0             ! vp in Pa
930 : 
931 : !       * Loop through every hour of day, where ik=hour
932 :         do ik = 1, 24
933 :           ii = in * 24 + ik - 24
934 : !         * (derived from 3.52+3.53) (Out[38], accounts for diurnal
935 : !           variation in air temperature
936 :           tair_h(ii) = tairmean + dtair * (0.0138d0                           &
937 :      &               * COS(3.513d0 - ((-1.d0 + ik) * p_pi) / 3.d0) + 0.0168d0 &
938 :      &               * COS(0.822d0 - ((-1.d0 + ik) * p_pi) / 4.d0) + 0.0984d0 &
939 :      &               * COS(0.360d0 - ((-1.d0 + ik) * p_pi) / 6.d0) + 0.4632d0 &
940 :      &               * COS(3.805d0 - ((-1.d0 + ik) * p_pi) / 12.d0))
941 : 
942 :           ca_h(ii) = ca_d(in) / 1.0d6
943 : 
944 : !         vd_h(ii) = 0.006028127d0 * 2.718282d0 ** ((17.27d0 * tair_h(ii)) &
945 : !    &             / (237.3d0 + tair_h(ii))) - 9.869233d-6 * vp__
946 : !         * (derived from 3.54+3.55) (Out[52]), accounts for diurnal variation in vapour deficit
947 :           vd_h(ii) = (((0.6108d0 * p_E ** (17.27d0 * tair_h(ii)        &
948 :      &             / (tair_h(ii) + 237.3d0))) * 1000) - vp__)          &
949 :      &             / (press_d(in) * 100.d0)
950 :           if (vd_h(ii) .le. 0.d0) vd_h(ii) = 0.d0
951 : 
952 : !         * average rainfall in hour ii (m/s)
953 :           rain_h(ii) = rain_d(in) / (24.d0 * 3600.d0 * 1000.d0)
954 : 
955 :           if (sunr .le. ik .and. ik + 1 .le. suns) then
956 : !           * ([Out30]), in mol/m2/s (derived from 3.51) accounts for
957 : !             diurnal variation in global irradiance
958 :             par_h(ii) = (-0.000873d0 * par_d(in) * COS(0.017453d0 * i_lat) &
959 :      &                * SQRT(0.920818d0 - 0.079182d0                   &
960 :      &                * COS(0.034428d0 * (10.d0 + dayyear(in))))       &
961 :      &                * COS(0.2618d0 * ik) - 0.000347d0 * par_d(in)    &
962 :      &                * COS(0.017214d0 * (10.d0 + dayyear(in)))        &
963 :      &                * SIN(0.017453d0 * i_lat)) / (-1.250192d0 * daylength &
964 :      &                * COS(0.017214d0 * (10.d0 + dayyear(in)))        &
965 :      &                * SIN(0.017453d0 * i_lat) + 24.d0                &
966 :      &                * COS(0.017453d0 * i_lat)                        &
967 :      &                * SQRT(0.920818d0 - 0.079182d0                   &
968 :      &                * COS(0.034428d0 * (10.d0 + dayyear(in))))       &
969 :      &                * (1.d0 - (0.158363d0                            &
970 :      &                * COS(0.017214d0 * (10.d0 + dayyear(in))) ** 2.d0 &
971 :      &                * TAN(0.017453d0 * i_lat) ** 2.d0)               &
972 :      &                / (0.920818d0 - 0.079182d0 * COS(0.034428d0      &
973 :      &                * (10.d0 + dayyear(in))))) ** 0.5d0)
974 :             if (par_h(ii) .lt. 0.d0) par_h(ii) = 0.d0
975 :           else
976 :             par_h(ii) = 0.d0
977 :           endif
978 : 
979 :           if (i_write_h == 1) then
980 :             write(kfile_hourlyweather,'(5i8,5e11.3)') ik, dayyear(in), &
981 :      &        fday(in), fmonth(in), fyear(in), tair_h(ii), vd_h(ii),   &
982 :      &        par_h(ii), rain_h(ii), ca_h(ii)
983 :           endif
984 : 
985 :         enddo
986 :       enddo
987 : 
988 :       if (i_write_h == 1) close(kfile_hourlyweather)
989 : 
990 :       return
991 :       end subroutine vom_calc_derived