!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! Calculates ref solution by integrating the model equation ! ! forward in time subject the boundary conditions bndcond and the ! ! forcing rhs. The routine is used to calculate the first guess, the ! ! representers, the reference case and the "final" solution, i.e., ! ! inverse estimate. ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! module m_reference contains subroutine reference(ref,bguess,MM,dz,depth,dt,tfin,vars,bnd0,bndH,chov) use mod_dimensions use mod_states use mod_variances use m_integrate use m_pseudo implicit none type(states), intent(out) :: ref(ndim) ! Model variable type(states), intent(out) :: bguess ! Model variable (init.values) type(variances), intent(in) :: vars ! Model variable type(local_states), intent(in) :: bnd0,bndH ! Boundary conditions real, intent(in) :: MM(ndim) ! Mixedlayer depth real, intent(in) :: dz ! Grid spacing in depth real, intent(in) :: depth ! Total depth real, intent(in) :: dt ! Grid spacing in time real, intent(in) :: tfin ! Final time real, intent(in) :: chov(kdim,kdim) ! Sqrt of spatial covariance matrix ! Automatic "arrays" type(states) forw ! "high res" sol real workA(kdim,antvar) real t0 ! Start time for integration real t1 ! End time for integration integer n,i,k,j character*2 tag ! initial cond. forw%N=10.0 forw%P=0.1 forw%H=0.1 ! Spinup period t0=0.0 ! t1=365.0 t1=1825.0 print '(a,2f10.2)','reference: Starting spinup',t0,t1 call integrate(forw,MM,bnd0,bndH,depth,t0,t1,dt) print '(a,2f10.2)','reference: Spinup ok',t0,t1 ! For the best guess: call pseudo(worka,kdim,antvar,chov) bguess%N=max(0.0,forw%N+sqrt(vars%ini%N)*worka(:,1)) bguess%P=max(0.0,forw%P+sqrt(vars%ini%P)*worka(:,2)) bguess%H=max(0.0,forw%H+sqrt(vars%ini%H)*worka(:,3)) ! For the reference solution: ! ref(1)=forw call pseudo(worka,kdim,antvar,chov) ref(1)%N=max(0.0,forw%N+sqrt(vars%ini%N)*worka(:,1)) ref(1)%P=max(0.0,forw%P+sqrt(vars%ini%P)*worka(:,2)) ref(1)%H=max(0.0,forw%H+sqrt(vars%ini%H)*worka(:,3)) do n=1,ndim-1 t0=float(n-1)*dt t1=float(n)*dt call integrate(forw,MM,bnd0,bndH,depth,t0,t1,dt) call pseudo(worka,kdim,antvar,chov) ref(n+1)%N=max(0.0, forw%N+sqrt(dt)*sqrt(vars%dyn%N)*worka(:,1)) ref(n+1)%P=max(0.0, forw%P+sqrt(dt)*sqrt(vars%dyn%P)*worka(:,2)) ref(n+1)%H=max(0.0, forw%H+sqrt(dt)*sqrt(vars%dyn%H)*worka(:,3)) ! ref(n+1)=forw enddo print '(a)','reference: Reference solution done' end subroutine reference end module