Program vom

  1 : !***********************************************************************
  2 : !        Optimised Vegetation Optimality Model (VOM)
  3 : !        Core program to run optimisation (sce) and transpmodel
  4 : !-----------------------------------------------------------------------
  5 : !        Author: Stan Schymanski, CWR, University of Western Australia
  6 : !        05/05/2004
  7 : !
  8 : !        Now at: MPI for Biogeochemistry, Jena, Germany
  9 : !        30/07/2007
 10 : !   sschym@bgc-jena.mpg.de
 11 : !
 12 : !-----------------------------------------------------------------------
 13 : !
 14 : !  Copyright (C) 2008  Stan Schymanski
 15 : !
 16 : !    This program is free software: you can redistribute it and/or modify
 17 : !    it under the terms of the GNU General Public License as published by
 18 : !    the Free Software Foundation, either version 3 of the License, or
 19 : !    (at your option) any later version.
 20 : !
 21 : !    This program is distributed in the hope that it will be useful,
 22 : !    but WITHOUT ANY WARRANTY; without even the implied warranty of
 23 : !    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 24 : !    GNU General Public License for more details.
 25 : !
 26 : !    You should have received a copy of the GNU General Public License
 27 : !    along with this program.  If not, see .
 28 : !
 29 : !***********************************************************************
 30 : 
 31 :       program vom
 32 :       use vom_sce_mod
 33 :       implicit none
 34 : 
 35 :       REAL*8, ALLOCATABLE :: vom_invar(:)
 36 :       REAL*8              :: vom_objfun
 37 : 
 38 :       INTEGER             :: beststat
 39 :       INTEGER             :: iostat
 40 :       LOGICAL             :: exist
 41 :       LOGICAL             :: run_sce
 42 :       REAL                :: starttime
 43 :       REAL                :: currtime
 44 :       REAL                :: runtime
 45 :       CHARACTER*100       :: outputpath_tmp ! Temporary outputpath 
 46 :       CHARACTER*100       :: inputpath_tmp  ! Temporary inputpath
 47 :       LOGICAL             :: change_in      ! Change input true/false
 48 :       LOGICAL             :: change_out     ! Change output true/false
 49 : 
 50 :       beststat = 0
 51 :       call cpu_time(starttime)
 52 : 
 53 :       call read_commandline(outputpath_tmp, inputpath_tmp, change_in, change_out)
 54 : 
 55 : !     * Parameter definitions
 56 : 
 57 :       call read_shufflepar()
 58 : 
 59 :       allocate(vom_invar(vom_npar))
 60 : 
 61 : 
 62 :       !optimize with sce
 63 :       if (vom_command .eq. 1 ) then
 64 : 
 65 :          call transpmodel_init_once(vom_command)
 66 : 
 67 :          if(sce_restart .eqv. .FALSE.) then 
 68 :             !remove old outputfiles for safety 
 69 :             open( kfile_sceout, FILE=trim(adjustl(i_outputpath)) // &
 70 :              trim(adjustl(sfile_sceout)), IOSTAT=iostat)
 71 :             if (iostat .eq. 0) close(kfile_sceout, status='delete')
 72 :             open( kfile_progress, FILE=trim(adjustl(i_outputpath)) // &
 73 :              trim(adjustl(sfile_progress)), STATUS='old', IOSTAT=iostat)
 74 :             if (iostat .eq. 0) close(kfile_progress, status='delete')
 75 :             open( kfile_lastloop, FILE=trim(adjustl(i_outputpath)) // &
 76 :              trim(adjustl(sfile_lastloop)), STATUS='old', IOSTAT=iostat)
 77 :             if (iostat .eq. 0) close(kfile_lastloop, status='delete')
 78 :             open( kfile_lastbest, FILE=trim(adjustl(i_outputpath)) // &
 79 :              trim(adjustl(sfile_lastbest)), STATUS='old', IOSTAT=iostat)
 80 :             if (iostat .eq. 0) close(kfile_lastbest, status='delete')
 81 :             open( kfile_bestpars, FILE=trim(adjustl(i_outputpath)) // &
 82 :              trim(adjustl(sfile_bestpars)), STATUS='old', IOSTAT=iostat)
 83 :             if (iostat .eq. 0) close(kfile_bestpars, status='delete')
 84 :             open( kfile_beststat, FILE=trim(adjustl(i_outputpath)) // &
 85 :              trim(adjustl(sfile_beststat)), STATUS='old', IOSTAT=iostat)
 86 :             if (iostat .eq. 0) close(kfile_beststat, status='delete')
 87 : 
 88 :             !run with initial seed
 89 :             write(*,*) "Initializing..."
 90 :             call sce_main()
 91 : 
 92 :          else
 93 :             write(*,*) "Restarting SCE from previous run..."
 94 : 
 95 :             !check if there are files to restart from
 96 :             open( kfile_lastloop, FILE=trim(adjustl(i_outputpath)) // &
 97 :              trim(adjustl(sfile_lastloop)), STATUS='old', IOSTAT=iostat)
 98 : 
 99 :             if (iostat .ne. 0) then
100 :                write(*,*) "ERROR: cannot restart, previous SCE-run did not finish a loop"
101 :                write(*,*) sfile_lastloop, " not found"
102 :                write(*,*) "Starting from scratch"
103 : 
104 :                !run with initial seed
105 :                write(*,*) "Initializing..."
106 :                call sce_main()
107 : 
108 :             end if
109 : 
110 :             !check if SCE already converged
111 :             open( kfile_beststat, FILE=trim(adjustl(i_outputpath)) // &
112 :              trim(adjustl(sfile_beststat)), STATUS='old', IOSTAT=iostat)
113 :             if (iostat .eq. 0) then
114 :                write(*,*) "ERROR: cannot restart from this state:"
115 :                write(*,*) "SCE already converged"
116 :                stop
117 :             end if
118 : 
119 :          end if 
120 : 
121 :          !run again untill sce_status.txt is written
122 :          run_sce = .TRUE.
123 :          do while (run_sce)
124 :             write(*,*) "Start looping"
125 :             call sce_main()
126 :             open(kfile_beststat, FILE=trim(adjustl(i_outputpath)) // &
127 :              trim(adjustl(sfile_beststat)), STATUS='old', IOSTAT=iostat)
128 :             
129 :             !stop if there is convergence  
130 :             if (iostat .eq. 0) then
131 :                run_sce = .FALSE.
132 :             end if
133 :             
134 :             call cpu_time(currtime) 
135 :             runtime = (currtime - starttime)  / (60.0) !minutes
136 :             !also stop when time is exceeded
137 :             if (runtime .gt. runtime_limit) then
138 :                run_sce = .FALSE.
139 :             end if
140 : 
141 :          end do
142 : 
143 :       endif
144 : 
145 :       !run normally
146 :       if (vom_command .eq. 2 ) then
147 : 
148 :          write(*,*) "Start single calculation with parameters..."
149 : 
150 :          call transpmodel_init_once(vom_command)
151 : 
152 :          open(kfile_pars, FILE=trim(adjustl(i_inputpath)) // &
153 :              trim(adjustl(sfile_pars)),&
154 :               STATUS='old', IOSTAT=iostat)
155 :               if (iostat .ne. 0) then
156 :                 write(0,*) "ERROR opening ", sfile_pars
157 :                 stop
158 :               endif
159 :               rewind(kfile_pars)
160 :               read(kfile_pars,*) vom_invar(:)
161 :             close(kfile_pars)
162 : 
163 : 
164 :         call transpmodel(vom_invar, SIZE(vom_invar), vom_objfun, vom_command)
165 : 
166 :       endif
167 : 
168 :       !run normally, no output except for ncp
169 :       if (vom_command .eq. 3 ) then
170 : 
171 :          write(*,*) "Start calculation of ncp with parameters..."
172 : 
173 :          open(kfile_pars, FILE=trim(adjustl(i_inputpath)) // &
174 :              trim(adjustl(sfile_pars)), STATUS='old', IOSTAT=iostat)
175 :               if (iostat .ne. 0) then
176 :                 write(0,*) "ERROR opening ", sfile_pars
177 :                 stop
178 :               endif
179 :               rewind(kfile_pars)
180 :               read(kfile_pars,*) vom_invar(:)
181 :             close(kfile_pars)
182 : 
183 :         call transpmodel_init_once(vom_command)
184 :         call transpmodel(vom_invar, SIZE(vom_invar), vom_objfun, vom_command)
185 : 
186 :         write(*,*) "The best carbon profit was: ",vom_objfun
187 :       endif
188 : 
189 : 
190 :       if (vom_command .eq. 5 ) then
191 :               write(*,*) "Random sampling of parameters ... "
192 :          call random_samples()
193 : 
194 :       endif
195 : 
196 : 
197 : 
198 :       deallocate(vom_invar)
199 : 
200 :       write(*,*) "Program terminated"
201 : 
202 :       end