/[MITgcm]/MITgcm_contrib/mlosch/optim_m1qn3/optim_sub.F
ViewVC logotype

Diff of /MITgcm_contrib/mlosch/optim_m1qn3/optim_sub.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.1 by mlosch, Thu Apr 26 11:10:06 2012 UTC revision 1.10 by mlosch, Tue Feb 12 15:50:49 2019 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4  c     Include ECCO_CPPOPTIONS because the ecco_ctrl,cost files  C     ECCO_CPPOPTIONS used to affect maxcvars and defined ALLOW_OBCS?_CONTROL
5  c     have headers with options for OBCS masks.  C#include "ECCO_CPPOPTIONS.h"
6  #include "ECCO_CPPOPTIONS.h"  C     now:
7    C     CTRL_OPTIONS affects maxcvars and may define ALLOW_OBCS?_CONTROL
8    #include "CTRL_OPTIONS.h"
9    
10        subroutine optim_sub(        subroutine optim_sub(
11       I                 nn, ff       I                 nn, ff
# Line 35  c     == global variables == Line 37  c     == global variables ==
37  #include "EEPARAMS.h"  #include "EEPARAMS.h"
38  #include "SIZE.h"  #include "SIZE.h"
39    
40    #if (defined (ALLOW_GENARR2D_CONTROL) || defined (ALLOW_GENARR3D_CONTROL) || defined (ALLOW_GENTIM2D_CONTROL))
41    # include "CTRL_SIZE.h"
42    #endif
43  #include "ctrl.h"  #include "ctrl.h"
44  #include "optim.h"  #include "optim.h"
45  #include "m1qn3_common.h"  #include "m1qn3_common.h"
# Line 48  c     == local variables == Line 53  c     == local variables ==
53    
54        _RL   objf        _RL   objf
55    
56  #if defined (DYNAMIC)  #ifdef DYNAMIC
57        _RL   xx(nn)        _RL, dimension(:), allocatable :: xx, adxx
       _RL   adxx(nn)  
 #elif defined (USE_POINTER) || (MAX_INDEPEND == 0)  
       _RL   xx  
       _RL   adxx  
       pointer (pxx,xx(1))  
       pointer (padxx,adxx(1))  
58  #else  #else
59        integer nmax        integer nmax
60        parameter( nmax = MAX_INDEPEND )        parameter( nmax = MAX_INDEPEND )
61        _RL   xx(nmax)        _RL   xx(nmax)
62        _RL   adxx(nmax)        _RL   adxx(nmax)
63  #endif  #endif
64          _RL  xxmean
65    
66    CML      logical coldStart
67  c     formal parameters of m1qn3  c     formal parameters of m1qn3
68        integer reverse        integer reverse
69        integer impres,imode(3),omode,niter,nsim,iz(5),indic        integer impres,imode(3),omode,niter,nsim,iz(5),indic
70        _RL dxmin,df1        _RL dxmin,df1
71        character*3 normtype        character*3 normtype
72  c     work arrays  c     work arrays
73        integer ndz, mupdate        integer ndz
74  CML      _RL dz(ndz)  CML      _RL dz(ndz)
75        double precision, dimension(:), allocatable :: dz        double precision, dimension(:), allocatable :: dz
76  c     extra dummy variables  c     extra dummy variables
# Line 90  c     == end of interface == Line 91  c     == end of interface ==
91    
92  c--   Allocate memory for the control variables and the gradient vector.  c--   Allocate memory for the control variables and the gradient vector.
93  #if defined(DYNAMIC)  #if defined(DYNAMIC)
94  #elif defined(USE_POINTER) || (MAX_INDEPEND == 0)        allocate(   xx(nn) )
95        call myalloc( pxx  ,  nn*REAL_BYTE )        allocate( adxx(nn) )
       call myalloc( padxx,  nn*REAL_BYTE )  
96  #endif  #endif
97    
98  #if defined (DYNAMIC)  #ifndef DYNAMIC
 #elif defined(USE_POINTER) || (MAX_INDEPEND == 0)  
 #else  
99        if (nn .gt. nmax) then        if (nn .gt. nmax) then
100          print*,' OPTIMUM: Not enough space.'          print*,' OPTIMUM: Not enough space.'
101          print*,'          nmax = ',nmax          print*,'          nmax = ',nmax
# Line 130  c Line 128  c
128        imode=(/0,1,0/)        imode=(/0,1,0/)
129        omode=-1        omode=-1
130  c     initialise work array  c     initialise work array
131        ndz = 3*nn+nupdate*(3*nn+1)        ndz = 4*nn+nupdate*(2*nn+1)
132        do i=1,5        do i=1,5
133         iz(i)=0         iz(i)=0
134        enddo        enddo
# Line 146  c     initialise the dummy arguments tha Line 144  c     initialise the dummy arguments tha
144        rzs(1)=UNSET_RS        rzs(1)=UNSET_RS
145        dzs(1)=UNSET_RL        dzs(1)=UNSET_RL
146    
147        if ( optimcycle .eq. 0 ) then  c--   first read the model output into xx, adxx, and cost function
148    c     value into objf
149          do i = 1,nn
150           xx(i)   = 0.
151           adxx(i) = 0.
152          enddo
153    c
154          print *, ' OPTIM_SUB: read model state'
155          call optim_readdata( nn, ctrlname, .false., objf,   xx )
156          call optim_readdata( nn, costname, .false., objf, adxx )
157          print *, ' OPTIM_SUB after reading ',
158         &           ctrlname, ' and ', costname, ':'
159          print *, ' OPTIM_SUB      nn = ', nn
160          print *, ' OPTIM_SUB    objf = ', objf
161          print *, ' OPTIM_SUB   xx(1) = ', xx(1)
162          print *, ' OPTIM_SUB adxx(1) = ', adxx(1)
163          
164          if ( coldStart ) then
165  c--   cold start  c--   cold start
166         print *, ' OPTIM_SUB: cold start, optimcycle =', optimcycle         print *, ' OPTIM_SUB: cold start, optimcycle =', optimcycle
167         imode(2) = 0         imode(2) = 0
# Line 154  c     this variable is the only one of t Line 169  c     this variable is the only one of t
169  c     that needs to be initialized here to make sure that we have a  c     that needs to be initialized here to make sure that we have a
170  c     clean start  c     clean start
171         reentry  = 0         reentry  = 0
172  c     ff has be read in optim_readparms, so we do not read it here again  c     compute expected decrease of cost function from objf and fmin;
173         objf = ff  c     this value is only used for a cold start of m1qn3_offline, for a
174         df1  = objf-fmin  c     warm start df1 is overwritten with data from a restart file
175           df1=objf-fmin
176           if ( df1 .le. 0. ) then
177            print *, ' OPTIM_SUB: df1 = objf-fmin = ', df1,
178         &       ' but should be > 0.'
179            stop 'ABNORMAL in S/R OPTIM_SUB'
180           endif
181    
182  c     open output file for m1qn3  c     open output file for m1qn3
183         open(io,file=fname_m1qn3,status='unknown')         open(io,file=fname_m1qn3,status='unknown')
184        else        else
# Line 170  c     requires restoring the state of m1 Line 192  c     requires restoring the state of m1
192  c     re-open output file for m1qn3  c     re-open output file for m1qn3
193         open(io,file=fname_m1qn3,status='old',position='append')         open(io,file=fname_m1qn3,status='old',position='append')
194        endif        endif
 c--   read the model output into xx,adxx  
       if ( indic .eq. 4 ) then  
        do i = 1,nn  
         xx(i)   = 0.  
         adxx(i) = 0.  
        enddo  
 c  
        print *, ' OPTIM_SUB: read model state'  
        call optim_readdata( nn, ctrlname, .false., objf,   xx )  
        call optim_readdata( nn, costname, .false., objf, adxx )  
        print *, ' OPTIM_SUB after reading nn, objf = ', nn, objf,  
      &      xx(1), adxx(1)  
       else  
        print *, ' OPTIM_SUB: indic = ', indic, ' is not possible'  
        stop 'ABNORMAL in S/R OPTIM_SUB'  
       endif  
195    
196  c--   call the minimizer, a slightly modified version of m1qn3 v3.3  c--   call the minimizer, a slightly modified version of m1qn3 v3.3
197  c     (Gilbert & Lemarechal, 1989), downloaded in April 2012.  c     (Gilbert & Lemarechal, 1989), downloaded in April 2012.
# Line 193  c     simul_rc,euclid,ctonbe,ctcabe are Line 199  c     simul_rc,euclid,ctonbe,ctcabe are
199  c     are provided within m1qn3. euclid, ctonbe, ctcabe can be replaced  c     are provided within m1qn3. euclid, ctonbe, ctcabe can be replaced
200  c     by something more efficient, simul_rc is a dummy routine for  c     by something more efficient, simul_rc is a dummy routine for
201  c     the reverse communication mode and should not be changed.  c     the reverse communication mode and should not be changed.
202        print *, ' OPTIM_SUB: call m1qn3_offline'        print *, ' OPTIM_SUB: call m1qn3_offline ........'
203        call m1qn3_offline (simul_rc,euclid,ctonbe,ctcabe,        call m1qn3_offline (simul_rc,euclid,ctonbe,ctcabe,
204       &     nn,xx,objf,adxx,dxmin,df1,       &     nn,xx,objf,adxx,dxmin,df1,
205       &     epsg,normtype,impres,io,imode,omode,niter,nsim,       &     epsg,normtype,impres,io,imode,omode,niter,nsim,
206       &     iz,dz,ndz,reverse,indic,izs,rzs,dzs)       &     iz,dz,ndz,reverse,indic,izs,rzs,dzs)
207        close(io)        close(io)
208          print *, ' OPTIM_SUB: ...........................'
209        print *, ' OPTIM_SUB: returned from m1qn3_offline'        print *, ' OPTIM_SUB: returned from m1qn3_offline'
210          print *, ' OPTIM_SUB:      nn = ', nn
211          print *, ' OPTIM_SUB:   xx(1) = ', xx(1), xx(2)
212          print *, ' OPTIM_SUB: adxx(1) = ', adxx(1), adxx(2)
213          print *, ' OPTIM_SUB: omode   = ', omode
214          print *, ' OPTIM_SUB: niter   = ', niter
215          print *, ' OPTIM_SUB: nsim    = ', nsim
216          print *, ' OPTIM_SUB: reverse = ', reverse
217    
218    c     compute min/max/mean/std of output vector see if it is within
219    c     reasonable bounds (prior to scaling)
220          xxmean = sum(xx)/dble(nn)
221          print *
222          print *,' OPTIM_SUB: mean(xx) =', xxmean
223          print *,' OPTIM_SUB:  max(xx) =', maxval(xx)
224          print *,' OPTIM_SUB:  min(xx) =', minval(xx)
225          print *,' OPTIM_SUB:  std(xx) =',sum((xx-xxmean)**2)/dble(nn)
226          print *
227  c     write state of m1qn3 into pickup file for warm restart  c     write state of m1qn3 into pickup file for warm restart
228        call optim_store_m1qn3(ndz,iz,dz,niter,nsim,epsg,df1,        call optim_store_m1qn3(ndz,iz,dz,niter,nsim,epsg,df1,
229       I     optimcycle,       I     optimcycle,
230       I     .true.)       I     .true.)
231  c     write model control vector  c     write model control vector
232        print *,' OPTIMS_SUB: writing ', nn,' sized control to file ',        print *,' OPTIM_SUB: writing ', nn,' sized control to file ',
233       &     ctrlname       &     ctrlname
234  c     give the cost function a funny value to make sure that nobody  c     give the cost function a funny value to make sure that nobody
235  c     mistakes it for the real one  c     mistakes it for the real one
236        call optim_writedata( nn, ctrlname, .false., -9999., xx )        call optim_writedata( nn, ctrlname, .false., -9999. _d 0, xx )
237    
238  c     clean up  c     clean up
239    #ifdef DYNAMIC
240          deallocate(xx, adxx)
241    #endif /* DYNAMIC */
242        deallocate(dz)        deallocate(dz)
243    
244  c     stopping criterion  c     stopping criterion
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.22