Some small bugs in feature/eCAI are fixed

.gitignore

@@ -24,8 +24,8 @@ bin/
 *.exe
 
 # others
-from_NR/*.f90
-downscale/numrec.f90
+#from_NR/*.f90
+#downscale/numrec.f90
 from_NR.tgz
 
 # no text files (OK)

downscale/Makefile-mac-tgq

@@ -0,0 +1,62 @@
+# ============= Makefile for downscaling program ======================
+
+# NOTE: first link one numrec.f file needed into this directory if not present
+#  eg, ln -s ../from_NR/numrec.f90
+
+# name of executable
+EXE = downscale.exe
+
+# core directory
+KOREDIR = /Users/localuser/Github/GMET/downscale
+
+# fortran compiler
+FC = gfortran
+
+# NetCDF libraries
+PATHNETCDF = /usr/local/Cellar/netcdf/4.6.3_1
+LIBNETCDF = -L$(PATHNETCDF)/lib -lnetcdff
+INCNETCDF = -I$(PATHNETCDF)/include
+
+# define compiler flags
+FLAGS = -p -g -Wall -ffree-line-length-none -fmax-errors=0 -fbacktrace -fcheck=bounds
+
+# define subroutines
+GMET_SUB = \
+  $(KOREDIR)/type.f90 \
+  $(KOREDIR)/precmod.f90 \
+  $(KOREDIR)/string_mod.f90 \
+  $(KOREDIR)/numrec.f90 \
+  $(KOREDIR)/utim.f90 \
+  $(KOREDIR)/read_config.f90 \
+  $(KOREDIR)/read_control.f90 \
+  $(KOREDIR)/read_nc.f90 \
+  $(KOREDIR)/read_domain_grid.f90 \
+  $(KOREDIR)/write_nc.f90 \
+    $(KOREDIR)/stats_routines.f90 \
+    $(KOREDIR)/regression_routines.f90 \
+    $(KOREDIR)/station_weights.f90 \
+    $(KOREDIR)/estimate_coefficients.f90 \
+    $(KOREDIR)/estimate_forcing_regression.f90 \
+    $(KOREDIR)/estimate_climo_regression.f90 \
+    $(KOREDIR)/estimate_climo_anom_regression.f90
+
+# Define the driver routine
+DRIVER = $(KOREDIR)/main.f90
+
+# Compile
+
+# tasks
+all: compile install clean
+
+# compile
+compile:
+	$(FC) $(FLAGS) $(GMET_SUB) $(DRIVER) \
+	$(LIBNETCDF) $(INCNETCDF) -o $(EXE)
+
+install:
+
+# Remove object files
+clean:
+	@rm -f *.o
+	@rm -f *.mod
+	@rm -f *__genmod.f90

downscale/Makefile-plato-tgq

@@ -0,0 +1,63 @@
+# ============= Makefile for downscaling program ======================
+
+# NOTE: first link one numrec.f file needed into this directory if not present
+#  eg, ln -s ../from_NR/numrec.f90
+
+# name of executable
+EXE = downscale.exe
+
+# core directory
+KOREDIR = /home/gut428/GMET/GMET_scripts/downscale
+
+# fortran compiler
+FC = gfortran
+
+# NetCDF libraries
+# PATHNETCDF = /usr/local
+PATHNETCDF = /cvmfs/soft.computecanada.ca/easybuild/software/2017/avx/Compiler/gcc5.4/netcdf-fortran/4.4.4
+LIBNETCDF = -L$(PATHNETCDF)/lib -lnetcdff
+INCNETCDF = -I$(PATHNETCDF)/include
+
+# define compiler flags
+FLAGS = -p -g -Wall -ffree-line-length-none -fmax-errors=0 -fbacktrace -fcheck=bounds
+
+# define subroutines
+GMET_SUB = \
+  $(KOREDIR)/type.f90 \
+  $(KOREDIR)/precmod.f90 \
+  $(KOREDIR)/string_mod.f90 \
+  $(KOREDIR)/numrec.f90 \
+  $(KOREDIR)/utim.f90 \
+  $(KOREDIR)/read_config.f90 \
+  $(KOREDIR)/read_control.f90 \
+  $(KOREDIR)/read_nc.f90 \
+  $(KOREDIR)/read_domain_grid.f90 \
+  $(KOREDIR)/write_nc.f90 \
+    $(KOREDIR)/stats_routines.f90 \
+    $(KOREDIR)/regression_routines.f90 \
+    $(KOREDIR)/station_weights.f90 \
+    $(KOREDIR)/estimate_coefficients.f90 \
+    $(KOREDIR)/estimate_forcing_regression.f90 \
+    $(KOREDIR)/estimate_climo_regression.f90 \
+    $(KOREDIR)/estimate_climo_anom_regression.f90
+
+# Define the driver routine
+DRIVER = $(KOREDIR)/main.f90
+
+# Compile
+
+# tasks
+all: compile install clean
+
+# compile
+compile:
+	$(FC) $(FLAGS) $(GMET_SUB) $(DRIVER) \
+	$(LIBNETCDF) $(INCNETCDF) -o $(EXE)
+
+install:
+
+# Remove object files
+clean:
+	@rm -f *.o
+	@rm -f *.mod
+	@rm -f *__genmod.f90

downscale/downscale.exe.dSYM/Contents/Info.plist

@@ -0,0 +1,20 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<!DOCTYPE plist PUBLIC "-//Apple Computer//DTD PLIST 1.0//EN" "http://www.apple.com/DTDs/PropertyList-1.0.dtd">
+<plist version="1.0">
+	<dict>
+		<key>CFBundleDevelopmentRegion</key>
+		<string>English</string>
+		<key>CFBundleIdentifier</key>
+		<string>com.apple.xcode.dsym.downscale.exe</string>
+		<key>CFBundleInfoDictionaryVersion</key>
+		<string>6.0</string>
+		<key>CFBundlePackageType</key>
+		<string>dSYM</string>
+		<key>CFBundleSignature</key>
+		<string>????</string>
+		<key>CFBundleShortVersionString</key>
+		<string>1.0</string>
+		<key>CFBundleVersion</key>
+		<string>1</string>
+	</dict>
+</plist>

downscale/estimate_climo_anom_regression.f90

@@ -317,7 +317,7 @@ end subroutine heapsort
   allocate(w_pcp_red_loocv(sta_limit-1,sta_limit-1))
   allocate(x_red_t_loocv(sta_limit-1,xsize))
   allocate(w_temp_red_loocv(sta_limit-1,sta_limit-1))
-  allocate(Y_tmean_red_loocv(sta_limit),Y_trange_red_loocv(sta_limit))
+  allocate(Y_tmean_red_loocv(sta_limit-1),Y_trange_red_loocv(sta_limit-1))
   allocate(tmp_weight_arr(sta_limit,sta_limit))
 
   ! initializations
@@ -467,7 +467,7 @@ end subroutine heapsort
     do i = 1, nstns, 1
       y (i) = prcp_data (i, t)
       y_tmean (i) = (tair_data(1, i, t)+tair_data(2, i, t)) / 2.0d0
-      y_trange (i) = abs (tair_data(2, i, t)-tair_data(1, i, t))
+      y_trange (i) = tair_data(2, i, t)-tair_data(1, i, t)
     end do
 
     ! do power transformation on precip vector (AWW: consider alternate transforms)
@@ -605,7 +605,7 @@ end subroutine heapsort
             twx_red = matmul (tx_red, w_pcp_red)
             call logistic_regression (x_red(:,1:3), y_red, twx_red, yp_red, b) ! AJN
 
-            if(-dot_product(Z(g,:),B) < 25.) then
+            if(-dot_product(Z(g,1:3),B) < 25.) then
               pop (g, t) = real (1.0/(1.0+exp(-dot_product(z(g, 1:3), b))), kind(sp))
             else
               POP(g,t) = 0.0
@@ -676,7 +676,7 @@ end subroutine heapsort
           twx_red = matmul (tx_red, w_temp_red)
           call least_squares (x_red_t(:,1:3), y_tmean_red, twx_red, b)
 
-          tmean (g, t) = real (dot_product(z(g, :), b), kind(sp))
+          tmean (g, t) = real (dot_product(z(g, 1:3), b), kind(sp))
 
           errsum = 0.0
           wgtsum = 0.0
@@ -711,7 +711,7 @@ end subroutine heapsort
           twx_red = matmul (tx_red, w_temp_red)
           call least_squares (x_red_t(:,1:3), y_trange_red, twx_red, b)
 
-          trange (g, t) = real (dot_product(z(g, :), b), kind(sp))
+          trange (g, t) = real (dot_product(z(g, 1:3), b), kind(sp))
 
           errsum = 0.0
           wgtsum = 0.0

downscale/estimate_climo_regression.f90

@@ -318,7 +318,7 @@ end subroutine heapsort
   allocate(w_pcp_red_loocv(sta_limit-1,sta_limit-1))
   allocate(x_red_t_loocv(sta_limit-1,xsize))
   allocate(w_temp_red_loocv(sta_limit-1,sta_limit-1))
-  allocate(Y_tmean_red_loocv(sta_limit),Y_trange_red_loocv(sta_limit))
+  allocate(Y_tmean_red_loocv(sta_limit-1),Y_trange_red_loocv(sta_limit-1))
   allocate(tmp_weight_arr(sta_limit,sta_limit))
 
   ! initializations

downscale/estimate_forcing_regression.f90

@@ -322,7 +322,7 @@ end subroutine heapsort
   allocate(w_pcp_red_loocv(sta_limit-1,sta_limit-1))
   allocate(x_red_t_loocv(sta_limit-1,xsize))
   allocate(w_temp_red_loocv(sta_limit-1,sta_limit-1))
-  allocate(Y_tmean_red_loocv(sta_limit),Y_trange_red_loocv(sta_limit))
+  allocate(Y_tmean_red_loocv(sta_limit-1),Y_trange_red_loocv(sta_limit-1))
   allocate(tmp_weight_arr(sta_limit,sta_limit))
 
 

downscale/numrec.f90

@@ -0,0 +1,137 @@
+! This collection of routines from num_rec was put together as part of the
+! downscaling code effort.  NCAR 2013
+
+
+subroutine nrerror (string)
+  character (len=*), intent (in) :: string
+  write (*,*) 'nrerror: ', string
+  stop 'program terminated by nrerror'
+end subroutine nrerror
+
+function assert_eq (n1, n2, n3, string)
+  character (len=*), intent (in) :: string
+  integer, intent (in) :: n1, n2, n3
+  integer :: assert_eq
+  if (n1 == n2 .and. n2 == n3) then
+    assert_eq = n1
+  else
+    write (*,*) 'nrerror: an assert_eq failed with this tag:', string
+    stop 'program terminated by assert_eq3'
+  end if
+end function assert_eq
+
+subroutine swap_d (a, b)
+  use type
+  real (dp), dimension (:), intent (inout) :: a, b
+  real (dp), dimension (size(a)) :: dum
+  dum = a
+  a = b
+  b = dum
+end subroutine swap_d
+
+function imaxloc (arr)
+  use type
+  real (dp), dimension (:), intent (in) :: arr
+  integer (i4b) :: imaxloc
+  integer (i4b), dimension (1) :: imax
+  imax = maxloc (arr(:))
+  imaxloc = imax (1)
+end function imaxloc
+
+function outerprod (a, b)
+  use type
+  real (dp), dimension (:), intent (in) :: a, b
+  real (dp), dimension (size(a), size(b)) :: outerprod
+  outerprod = spread (a, dim=2, ncopies=size(b)) * spread (b, dim=1, ncopies=size(a))
+end function outerprod
+
+subroutine ludcmp (a, indx, d)
+  use type
+  implicit none
+  interface
+    subroutine nrerror (string)
+      character (len=*), intent (in) :: string
+    end subroutine nrerror
+
+    function assert_eq (n1, n2, n3, string)
+      character (len=*), intent (in) :: string
+      integer, intent (in) :: n1, n2, n3
+      integer :: assert_eq
+    end function assert_eq
+
+    function imaxloc (arr)
+      use type
+      real (dp), dimension (:), intent (in) :: arr
+      integer (i4b) :: imaxloc
+    end function imaxloc
+
+    function outerprod (a, b)
+      use type
+      real (dp), dimension (:), intent (in) :: a, b
+      real (dp), dimension (size(a), size(b)) :: outerprod
+    end function outerprod
+
+    subroutine swap_d (a, b)
+      use type
+      real (dp), dimension (:), intent (inout) :: a, b
+    end subroutine swap_d
+  end interface
+
+  real (dp), dimension (:, :), intent (inout) :: a
+  integer (i4b), dimension (:), intent (out) :: indx
+  real (dp), intent (out) :: d
+  real (dp), dimension (size(a, 1)) :: vv
+!  REAL(DP), PARAMETER :: TINY=1.0e-20
+  integer (i4b) :: j, n, imax
+  n = assert_eq (size(a, 1), size(a, 2), size(indx), 'ludcmp')
+  d = 1.0d0
+  vv = maxval (abs(a), dim=2)
+  if (any(vv == 0.0)) call nrerror ('singular matrix in ludcmp')
+  vv = 1.0d0 / vv
+  do j = 1, n
+    imax = (j-1) + imaxloc (vv(j:n)*abs(a(j:n, j)))
+    if (j /= imax) then
+      call swap_d (a(imax, :), a(j, :))
+      d = - d
+      vv (imax) = vv (j)
+    end if
+    indx (j) = imax
+    if (a(j, j) == 0.0) a (j, j) = tiny
+    a (j+1:n, j) = a (j+1:n, j) / a (j, j)
+    a (j+1:n, j+1:n) = a (j+1:n, j+1:n) - outerprod (a(j+1:n, j), a(j, j+1:n))
+  end do
+end subroutine ludcmp
+
+subroutine lubksb (a, indx, b)
+  use type
+  implicit none
+  interface
+    function assert_eq (n1, n2, n3, string)
+      character (len=*), intent (in) :: string
+      integer, intent (in) :: n1, n2, n3
+      integer :: assert_eq
+    end function assert_eq
+  end interface
+
+  real (dp), dimension (:, :), intent (in) :: a
+  integer (i4b), dimension (:), intent (in) :: indx
+  real (dp), dimension (:), intent (inout) :: b
+  integer (i4b) :: i, n, ii, ll
+  real (dp) :: summ
+  n = assert_eq (size(a, 1), size(a, 2), size(indx), 'lubksb')
+  ii = 0
+  do i = 1, n
+    ll = indx (i)
+    summ = b (ll)
+    b (ll) = b (i)
+    if (ii /= 0) then
+      summ = summ - dot_product (a(i, ii:i-1), b(ii:i-1))
+    else if (summ /= 0.0) then
+      ii = i
+    end if
+    b (i) = summ
+  end do
+  do i = n, 1, - 1
+    b (i) = (b(i)-dot_product(a(i, i+1:n), b(i+1:n))) / a (i, i)
+  end do
+end subroutine lubksb

from_NR/erf.f90

@@ -0,0 +1,20 @@
+	FUNCTION erf_s(x)
+	USE nrtype
+	USE nr, ONLY : gammp
+	IMPLICIT NONE
+	REAL(SP), INTENT(IN) :: x
+	REAL(SP) :: erf_s
+	erf_s=gammp(0.5_sp,x**2)
+	if (x < 0.0) erf_s=-erf_s
+	END FUNCTION erf_s
+
+
+	FUNCTION erf_v(x)
+	USE nrtype
+	USE nr, ONLY : gammp
+	IMPLICIT NONE
+	REAL(SP), DIMENSION(:), INTENT(IN) :: x
+	REAL(SP), DIMENSION(size(x)) :: erf_v
+	erf_v=gammp(spread(0.5_sp,1,size(x)),x**2)
+	where (x < 0.0) erf_v=-erf_v
+	END FUNCTION erf_v