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