m_reference.F90 2.92 KB
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                     !
!  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