308
  309
  310
  311
  312      integer, intent(in) :: ng, tile, model
  313      integer, intent(in) :: LBi, UBi, LBj, UBj, LBij, UBij
  314      integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
  315      integer, intent(in) :: Lold, Lnew, outLoop
  316
  317# ifdef ASSUMED_SHAPE
  318#  ifdef MASKING
  319      real(r8), intent(in) :: rmask(LBi:,LBj:)
  320      real(r8), intent(in) :: umask(LBi:,LBj:)
  321      real(r8), intent(in) :: vmask(LBi:,LBj:)
  322#  endif
  323#   ifdef ADJUST_BOUNDARY
  324#    ifdef SOLVE3D
  325      real(r8), intent(in) :: t_obc_std(LBij:,:,:,:)
  326      real(r8), intent(in) :: u_obc_std(LBij:,:,:)
  327      real(r8), intent(in) :: v_obc_std(LBij:,:,:)
  328#    endif
  329      real(r8), intent(in) :: ubar_obc_std(LBij:,:)
  330      real(r8), intent(in) :: vbar_obc_std(LBij:,:)
  331      real(r8), intent(in) :: zeta_obc_std(LBij:,:)
  332#   endif
  333#   ifdef ADJUST_WSTRESS
  334      real(r8), intent(in) :: sustr_std(LBi:,LBj:)
  335      real(r8), intent(in) :: svstr_std(LBi:,LBj:)
  336#   endif
  337#   if defined ADJUST_STFLUX && defined SOLVE3D
  338      real(r8), intent(in) :: stflx_std(LBi:,LBj:,:)
  339#   endif
  340#   ifdef SOLVE3D
  341      real(r8), intent(in) :: t_std(LBi:,LBj:,:,:,:)
  342      real(r8), intent(in) :: u_std(LBi:,LBj:,:,:)
  343      real(r8), intent(in) :: v_std(LBi:,LBj:,:,:)
  344#    if defined WEAK_CONSTRAINT && defined TIME_CONV
  345      real(r8), intent(in) :: ubar_std(LBi:,LBj:,:)
  346      real(r8), intent(in) :: vbar_std(LBi:,LBj:,:)
  347#    endif
  348#   else
  349      real(r8), intent(in) :: ubar_std(LBi:,LBj:,:)
  350      real(r8), intent(in) :: vbar_std(LBi:,LBj:,:)
  351#   endif
  352      real(r8), intent(in) :: zeta_std(LBi:,LBj:,:)
  353#  ifdef ADJUST_BOUNDARY
  354#   ifdef SOLVE3D
  355      real(r8), intent(inout) :: ad_t_obc(LBij:,:,:,:,:,:)
  356      real(r8), intent(inout) :: ad_u_obc(LBij:,:,:,:,:)
  357      real(r8), intent(inout) :: ad_v_obc(LBij:,:,:,:,:)
  358#   endif
  359      real(r8), intent(inout) :: ad_ubar_obc(LBij:,:,:,:)
  360      real(r8), intent(inout) :: ad_vbar_obc(LBij:,:,:,:)
  361      real(r8), intent(inout) :: ad_zeta_obc(LBij:,:,:,:)
  362#  endif
  363#  ifdef ADJUST_WSTRESS
  364      real(r8), intent(inout) :: ad_ustr(LBi:,LBj:,:,:)
  365      real(r8), intent(inout) :: ad_vstr(LBi:,LBj:,:,:)
  366#  endif
  367#  ifdef SOLVE3D
  368#   ifdef ADJUST_STFLUX
  369      real(r8), intent(inout) :: ad_tflux(LBi:,LBj:,:,:,:)
  370#   endif
  371      real(r8), intent(inout) :: ad_t(LBi:,LBj:,:,:,:)
  372      real(r8), intent(inout) :: ad_u(LBi:,LBj:,:,:)
  373      real(r8), intent(inout) :: ad_v(LBi:,LBj:,:,:)
  374#   if defined WEAK_CONSTRAINT && defined TIME_CONV
  375      real(r8), intent(inout) :: ad_ubar(LBi:,LBj:,:)
  376      real(r8), intent(inout) :: ad_vbar(LBi:,LBj:,:)
  377#   endif
  378#  else
  379      real(r8), intent(inout) :: ad_ubar(LBi:,LBj:,:)
  380      real(r8), intent(inout) :: ad_vbar(LBi:,LBj:,:)
  381#  endif
  382      real(r8), intent(inout) :: ad_zeta(LBi:,LBj:,:)
  383#  ifdef ADJUST_BOUNDARY
  384#   ifdef SOLVE3D
  385      real(r8), intent(inout) :: nl_t_obc(LBij:,:,:,:,:,:)
  386      real(r8), intent(inout) :: nl_u_obc(LBij:,:,:,:,:)
  387      real(r8), intent(inout) :: nl_v_obc(LBij:,:,:,:,:)
  388#   endif
  389      real(r8), intent(inout) :: nl_ubar_obc(LBij:,:,:,:)
  390      real(r8), intent(inout) :: nl_vbar_obc(LBij:,:,:,:)
  391      real(r8), intent(inout) :: nl_zeta_obc(LBij:,:,:,:)
  392#  endif
  393#  ifdef ADJUST_WSTRESS
  394      real(r8), intent(inout) :: nl_ustr(LBi:,LBj:,:,:)
  395      real(r8), intent(inout) :: nl_vstr(LBi:,LBj:,:,:)
  396#  endif
  397#  ifdef SOLVE3D
  398#   ifdef ADJUST_STFLUX
  399      real(r8), intent(inout) :: nl_tflux(LBi:,LBj:,:,:,:)
  400#   endif
  401      real(r8), intent(inout) :: nl_t(LBi:,LBj:,:,:,:)
  402      real(r8), intent(inout) :: nl_u(LBi:,LBj:,:,:)
  403      real(r8), intent(inout) :: nl_v(LBi:,LBj:,:,:)
  404#   if defined WEAK_CONSTRAINT && defined TIME_CONV
  405      real(r8), intent(inout) :: nl_ubar(LBi:,LBj:,:)
  406      real(r8), intent(inout) :: nl_vbar(LBi:,LBj:,:)
  407#   endif
  408#  else
  409      real(r8), intent(inout) :: nl_ubar(LBi:,LBj:,:)
  410      real(r8), intent(inout) :: nl_vbar(LBi:,LBj:,:)
  411#  endif
  412      real(r8), intent(inout) :: nl_zeta(LBi:,LBj:,:)
  413#  ifdef ADJUST_BOUNDARY
  414#   ifdef SOLVE3D
  415      real(r8), intent(inout) :: tl_t_obc(LBij:,:,:,:,:,:)
  416      real(r8), intent(inout) :: tl_u_obc(LBij:,:,:,:,:)
  417      real(r8), intent(inout) :: tl_v_obc(LBij:,:,:,:,:)
  418#   endif
  419      real(r8), intent(inout) :: tl_ubar_obc(LBij:,:,:,:)
  420      real(r8), intent(inout) :: tl_vbar_obc(LBij:,:,:,:)
  421      real(r8), intent(inout) :: tl_zeta_obc(LBij:,:,:,:)
  422#  endif
  423#  ifdef ADJUST_WSTRESS
  424      real(r8), intent(inout) :: tl_ustr(LBi:,LBj:,:,:)
  425      real(r8), intent(inout) :: tl_vstr(LBi:,LBj:,:,:)
  426#  endif
  427#  ifdef SOLVE3D
  428#   ifdef ADJUST_STFLUX
  429      real(r8), intent(inout) :: tl_tflux(LBi:,LBj:,:,:,:)
  430#   endif
  431      real(r8), intent(inout) :: tl_t(LBi:,LBj:,:,:,:)
  432      real(r8), intent(inout) :: tl_u(LBi:,LBj:,:,:)
  433      real(r8), intent(inout) :: tl_v(LBi:,LBj:,:,:)
  434#   if defined WEAK_CONSTRAINT && defined TIME_CONV
  435      real(r8), intent(inout) :: tl_ubar(LBi:,LBj:,:)
  436      real(r8), intent(inout) :: tl_vbar(LBi:,LBj:,:)
  437#   endif
  438#  else
  439      real(r8), intent(inout) :: tl_ubar(LBi:,LBj:,:)
  440      real(r8), intent(inout) :: tl_vbar(LBi:,LBj:,:)
  441#  endif
  442      real(r8), intent(inout) :: tl_zeta(LBi:,LBj:,:)
  443 
  444# else
  445 
  446#  ifdef MASKING
  447      real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
  448      real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
  449      real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
  450#  endif
  451#  ifdef ADJUST_BOUNDARY
  452#   ifdef SOLVE3D
  453      real(r8), intent(in) :: t_obc_std(LBij:UBij,N(ng),4,NT(ng))
  454      real(r8), intent(in) :: u_obc_std(LBij:UBij,N(ng),4)
  455      real(r8), intent(in) :: v_obc_std(LBij:UBij,N(ng),4)
  456#   endif
  457      real(r8), intent(in) :: ubar_obc_std(LBij:UBij,4)
  458      real(r8), intent(in) :: vbar_obc_std(LBij:UBij,4)
  459      real(r8), intent(in) :: zeta_obc_std(LBij:UBij,4)
  460#  endif
  461#  ifdef ADJUST_WSTRESS
  462      real(r8), intent(in) :: sustr_std(LBi:,LBj:)
  463      real(r8), intent(in) :: svstr_std(LBi:,LBj:)
  464#  endif
  465#  if defined ADJUST_STFLUX && defined SOLVE3D
  466      real(r8), intent(in) :: stflx_std(LBi:UBi,LBj:UBj,NT(ng))
  467#  endif
  468#  ifdef SOLVE3D
  469      real(r8), intent(in) :: t_std(LBi:UBi,LBj:UBj,N(ng),NSA,NT(ng))
  470      real(r8), intent(in) :: u_std(LBi:UBi,LBj:UBj,N(ng),NSA)
  471      real(r8), intent(in) :: v_std(LBi:UBi,LBj:UBj,N(ng),NSA)
  472#   if defined WEAK_CONSTRAINT && defined TIME_CONV
  473      real(r8), intent(in) :: ubar_std(LBi:UBi,LBj:UBj,NSA)
  474      real(r8), intent(in) :: vbar_std(LBi:UBi,LBj:UBj,NSA)
  475#   endif
  476#  else
  477      real(r8), intent(in) :: ubar_std(LBi:UBi,LBj:UBj,NSA)
  478      real(r8), intent(in) :: vbar_std(LBi:UBi,LBj:UBj,NSA)
  479#  endif
  480      real(r8), intent(in) :: zeta_std(LBi:UBi,LBj:UBj,NSA)
  481#  ifdef ADJUST_BOUNDARY
  482#   ifdef SOLVE3D
  483      real(r8), intent(inout) :: ad_t_obc(LBij:UBij,N(ng),4,            &
  484     &                                    Nbrec(ng),2,NT(ng))
  485      real(r8), intent(inout) :: ad_u_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
  486      real(r8), intent(inout) :: ad_v_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
  487#   endif
  488      real(r8), intent(inout) :: ad_ubar_obc(LBij:UBij,4,Nbrec(ng),2)
  489      real(r8), intent(inout) :: ad_vbar_obc(LBij:UBij,4,Nbrec(ng),2)
  490      real(r8), intent(inout) :: ad_zeta_obc(LBij:UBij,4,Nbrec(ng),2)
  491#  endif
  492#  ifdef ADJUST_WSTRESS
  493      real(r8), intent(inout) :: ad_ustr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
  494      real(r8), intent(inout) :: ad_vstr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
  495#  endif
  496#  ifdef SOLVE3D
  497#   ifdef ADJUST_STFLUX
  498      real(r8), intent(inout) :: ad_tflux(LBi:UBi,LBj:UBj,              &
  499     &                                    Nfrec(ng),2,NT(ng))
  500#   endif
  501      real(r8), intent(inout) :: ad_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
  502      real(r8), intent(inout) :: ad_u(LBi:UBi,LBj:UBj,N(ng),2)
  503      real(r8), intent(inout) :: ad_v(LBi:UBi,LBj:UBj,N(ng),2)
  504#   if defined WEAK_CONSTRAINT && defined TIME_CONV
  505      real(r8), intent(inout) :: ad_ubar(LBi:UBi,LBj:UBj,:)
  506      real(r8), intent(inout) :: ad_vbar(LBi:UBi,LBj:UBj,:)
  507#   endif
  508#  else
  509      real(r8), intent(inout) :: ad_ubar(LBi:UBi,LBj:UBj,:)
  510      real(r8), intent(inout) :: ad_vbar(LBi:UBi,LBj:UBj,:)
  511#  endif
  512      real(r8), intent(inout) :: ad_zeta(LBi:UBi,LBj:UBj,:)
  513#  ifdef ADJUST_BOUNDARY
  514#   ifdef SOLVE3D
  515      real(r8), intent(inout) :: nl_t_obc(LBij:UBij,N(ng),4,            &
  516     &                                    Nbrec(ng),2,NT(ng))
  517      real(r8), intent(inout) :: nl_u_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
  518      real(r8), intent(inout) :: nl_v_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
  519#   endif
  520      real(r8), intent(inout) :: nl_ubar_obc(LBij:UBij,4,Nbrec(ng),2)
  521      real(r8), intent(inout) :: nl_vbar_obc(LBij:UBij,4,Nbrec(ng),2)
  522      real(r8), intent(inout) :: nl_zeta_obc(LBij:UBij,4,Nbrec(ng),2)
  523#  endif
  524#  ifdef ADJUST_WSTRESS
  525      real(r8), intent(inout) :: nl_ustr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
  526      real(r8), intent(inout) :: nl_vstr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
  527#  endif
  528#  ifdef SOLVE3D
  529#   ifdef ADJUST_STFLUX
  530      real(r8), intent(inout) :: nl_tflux(LBi:UBi,LBj:UBj,              &
  531     &                                    Nfrec(ng),2,NT(ng))
  532#   endif
  533      real(r8), intent(inout) :: nl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
  534      real(r8), intent(inout) :: nl_u(LBi:UBi,LBj:UBj,N(ng),2)
  535      real(r8), intent(inout) :: nl_v(LBi:UBi,LBj:UBj,N(ng),2)
  536#   if defined WEAK_CONSTRAINT && defined TIME_CONV
  537      real(r8), intent(inout) :: nl_ubar(LBi:UBi,LBj:UBj,:)
  538      real(r8), intent(inout) :: nl_vbar(LBi:UBi,LBj:UBj,:)
  539#   endif
  540#  else
  541      real(r8), intent(inout) :: nl_ubar(LBi:UBi,LBj:UBj,:)
  542      real(r8), intent(inout) :: nl_vbar(LBi:UBi,LBj:UBj,:)
  543#  endif
  544      real(r8), intent(inout) :: nl_zeta(LBi:UBi,LBj:UBj,:)
  545#  ifdef ADJUST_BOUNDARY
  546#   ifdef SOLVE3D
  547      real(r8), intent(inout) :: tl_t_obc(LBij:UBij,N(ng),4,            &
  548     &                                    Nbrec(ng),2,NT(ng))
  549      real(r8), intent(inout) :: tl_u_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
  550      real(r8), intent(inout) :: tl_v_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
  551#   endif
  552      real(r8), intent(inout) :: tl_ubar_obc(LBij:UBij,4,Nbrec(ng),2)
  553      real(r8), intent(inout) :: tl_vbar_obc(LBij:UBij,4,Nbrec(ng),2)
  554      real(r8), intent(inout) :: tl_zeta_obc(LBij:UBij,4,Nbrec(ng),2)
  555#  endif
  556#  ifdef ADJUST_WSTRESS
  557      real(r8), intent(inout) :: tl_ustr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
  558      real(r8), intent(inout) :: tl_vstr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
  559#  endif
  560#  ifdef SOLVE3D
  561#   ifdef ADJUST_STFLUX
  562      real(r8), intent(inout) :: tl_tflux(LBi:UBi,LBj:UBj,              &
  563     &                                    Nfrec(ng),2,NT(ng))
  564#   endif
  565      real(r8), intent(inout) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
  566      real(r8), intent(inout) :: tl_u(LBi:UBi,LBj:UBj,N(ng),2)
  567      real(r8), intent(inout) :: tl_v(LBi:UBi,LBj:UBj,N(ng),2)
  568#   if defined WEAK_CONSTRAINT && defined TIME_CONV
  569      real(r8), intent(inout) :: tl_ubar(LBi:UBi,LBj:UBj,:)
  570      real(r8), intent(inout) :: tl_vbar(LBi:UBi,LBj:UBj,:)
  571#   endif
  572#  else
  573      real(r8), intent(inout) :: tl_ubar(LBi:UBi,LBj:UBj,:)
  574      real(r8), intent(inout) :: tl_vbar(LBi:UBi,LBj:UBj,:)
  575#  endif
  576      real(r8), intent(inout) :: tl_zeta(LBi:UBi,LBj:UBj,:)
  577# endif
  578
  579
  580
  581      integer :: L1 = 1
  582      integer :: L2 = 2
  583 
  584      integer :: Linp, Lout, Lscale, Lwrk, Lwrk1, i, j, ic
  585      integer :: info, ivec, jvec, rec
  586 
  587      real(r8) :: fac, fac1, fac2
  588 
  589      logical :: Lweak
  590 
  591      character (len=256) :: ncname
  592 
  593# include "set_bounds.h"
  594
  595      calledfrom=myfile
  596      sourcefile=myfile
  597
  598
  599
  600# if defined POSTERIOR_ERROR_I
  601
  602# elif defined POSTERIOR_ERROR_F
  603
  604# endif
  605
  606
  607
  608
  609
  610
  611
  612
  613      IF (master) WRITE (stdout,10)
  614 10   FORMAT (/,' <<<< Full Posterior Error Covariance Matrix >>>>',/)
  615
  616      lweak=.false.
  617
  618
  619
  620
  621      zlanczos_coef=0.0_r8
  622      DO i=1,ninner
  623        zlanczos_coef(i,i)=cg_delta(i,outloop)
  624      END DO
  625      DO i=1,ninner-1
  626        zlanczos_coef(i,i+1)=cg_beta(i+1,outloop)
  627      END DO
  628      DO i=2,ninner
  629        zlanczos_coef(i,i-1)=cg_beta(i,outloop)
  630      END DO
  631      zlanczos_inv=zlanczos_coef
  632
  633
  634
  635
  636
  637      IF (master) THEN
  638        CALL dpotrf ('U', ninner, zlanczos_inv, ninner, info)
  639      END IF
  640# ifdef DISTRIBUTE
  641      CALL mp_bcasti (ng, model, info)
  642# endif
  643      IF (info.ne.0) THEN
  644        IF (master) WRITE (stdout,*) ' Error in DPOTRF: info = ', info
  645        exit_flag=8
  646        RETURN
  647      END IF
  648
  649      IF (master) THEN
  650        CALL dpotri ('U', ninner, zlanczos_inv, ninner, info)
  651      END IF
  652# ifdef DISTRIBUTE
  653      CALL mp_bcasti (ng, model, info)
  654      CALL mp_bcastf (ng, model, zlanczos_inv)
  655# endif
  656      IF (info.ne.0) THEN
  657        IF (master) WRITE (stdout,*) ' Error in DPOTRI: info = ', info
  658        exit_flag=8
  659        RETURN
  660      END IF
  661
  662
  663
  664
  665      DO j=1,ninner
  666        DO i=j+1,ninner
  667          zlanczos_inv(i,j)=zlanczos_inv(j,i)
  668        END DO
  669      END DO
  670
  671
  672
  673      DO j=1,ninner
  674       DO i=1,ninner
  675         zlanczos_err(i,j)=0.0_r8
  676         DO ic=1,ninner
  677            zlanczos_err(i,j)=zlanczos_err(i,j)+                        &
  678     &                        zlanczos_inv(i,ic)*zlanczos_coef(ic,j)
  679         END DO
  680       END DO
  681      END DO
  682
  683
  684
  685
  686      fac=0.0_r8
  687      CALL state_initialize (ng, tile,                                  &
  688     &                       lbi, ubi, lbj, ubj, lbij, ubij,            &
  689     &                       l1, fac,                                   &
  690# ifdef MASKING
  691     &                       rmask, umask, vmask,                       &
  692# endif
  693# ifdef ADJUST_BOUNDARY
  694#  ifdef SOLVE3D
  695     &                       tl_t_obc, tl_u_obc, tl_v_obc,              &
  696#  endif
  697     &                       tl_ubar_obc, tl_vbar_obc,                  &
  698     &                       tl_zeta_obc,                               &
  699# endif
  700# ifdef ADJUST_WSTRESS
  701     &                       tl_ustr, tl_vstr,                          &
  702# endif
  703# ifdef SOLVE3D
  704#  ifdef ADJUST_STFLUX
  705     &                       tl_tflux,                                  &
  706#  endif
  707     &                       tl_t, tl_u, tl_v,                          &
  708# else
  709     &                       tl_ubar, tl_vbar,                          &
  710# endif
  711     &                       tl_zeta)
  712
  713
  714
  715
  716
  717      ncname=hss(ng)%name
  718      DO ivec=1,ninner
  719        rec=ivec
  720
  721
  722
  723
  724        CALL state_read (ng, tile, model, hss(ng)%IOtype,               &
  725     &                   lbi, ubi, lbj, ubj, lbij, ubij,                &
  726     &                   l1, rec,                                       &
  727     &                   0, hss(ng)%ncid, ncname,                       &
  728# ifdef MASKING
  729     &                   rmask, umask, vmask,                           &
  730# endif
  731# ifdef ADJUST_BOUNDARY
  732#  ifdef SOLVE3D
  733     &                   ad_t_obc, ad_u_obc, ad_v_obc,                  &
  734#  endif
  735     &                   ad_ubar_obc, ad_vbar_obc,                      &
  736     &                   ad_zeta_obc,                                   &
  737# endif
  738# ifdef ADJUST_WSTRESS
  739     &                   ad_ustr, ad_vstr,                              &
  740# endif
  741# ifdef SOLVE3D
  742#  ifdef ADJUST_STFLUX
  743     &                   ad_tflux,                                      &
  744#  endif
  745     &                   ad_t, ad_u, ad_v,                              &
  746# else
  747     &                   ad_ubar, ad_vbar,                              &
  748# endif
  749     &                   ad_zeta)
  750
  751
  752
  753
  754        CALL state_product (ng, tile, model,                            &
  755     &                      lbi, ubi, lbj, ubj, lbij, ubij,             &
  756# ifdef MASKING
  757     &                      rmask, umask, vmask,                        &
  758# endif
  759# ifdef ADJUST_BOUNDARY
  760#  ifdef SOLVE3D
  761     &                      ad_t_obc(:,:,:,:,l1,:),                     &
  762     &                      ad_t_obc(:,:,:,:,l1,:),                     &
  763     &                      nl_t_obc(:,:,:,:,l1,:),                     &
  764     &                      ad_u_obc(:,:,:,:,l1),                       &
  765     &                      ad_u_obc(:,:,:,:,l1),                       &
  766     &                      nl_u_obc(:,:,:,:,l1),                       &
  767     &                      ad_v_obc(:,:,:,:,l1),                       &
  768     &                      ad_v_obc(:,:,:,:,l1),                       &
  769     &                      nl_v_obc(:,:,:,:,l1),                       &
  770#  endif
  771     &                      ad_ubar_obc(:,:,:,l1),                      &
  772     &                      ad_ubar_obc(:,:,:,l1),                      &
  773     &                      nl_ubar_obc(:,:,:,l1),                      &
  774     &                      ad_vbar_obc(:,:,:,l1),                      &
  775     &                      ad_vbar_obc(:,:,:,l1),                      &
  776     &                      nl_vbar_obc(:,:,:,l1),                      &
  777     &                      ad_zeta_obc(:,:,:,l1),                      &
  778     &                      ad_zeta_obc(:,:,:,l1),                      &
  779     &                      nl_zeta_obc(:,:,:,l1),                      &
  780# endif
  781# ifdef ADJUST_WSTRESS
  782     &                      ad_ustr(:,:,:,l1),                          &
  783     &                      ad_ustr(:,:,:,l1),                          &
  784     &                      nl_ustr(:,:,:,l1),                          &
  785     &                      ad_vstr(:,:,:,l1),                          &
  786     &                      ad_vstr(:,:,:,l1),                          &
  787     &                      nl_vstr(:,:,:,l1),                          &
  788# endif
  789# ifdef SOLVE3D
  790#  ifdef ADJUST_STFLUX
  791     &                      ad_tflux(:,:,:,l1,:),                       &
  792     &                      ad_tflux(:,:,:,l1,:),                       &
  793     &                      nl_tflux(:,:,:,l1,:),                       &
  794#  endif
  795     &                      ad_t(:,:,:,l1,:),                           &
  796     &                      ad_t(:,:,:,l1,:),                           &
  797     &                      nl_t(:,:,:,l1,:),                           &
  798     &                      ad_u(:,:,:,l1),                             &
  799     &                      ad_u(:,:,:,l1),                             &
  800     &                      nl_u(:,:,:,l1),                             &
  801     &                      ad_v(:,:,:,l1),                             &
  802     &                      ad_v(:,:,:,l1),                             &
  803     &                      nl_v(:,:,:,l1),                             &
  804# else
  805     &                      ad_ubar(:,:,l1),                            &
  806     &                      ad_ubar(:,:,l1),                            &
  807     &                      nl_ubar(:,:,l1),                            &
  808     &                      ad_vbar(:,:,l1),                            &
  809     &                      ad_vbar(:,:,l1),                            &
  810     &                      nl_vbar(:,:,l1),                            &
  811# endif
  812     &                      ad_zeta(:,:,l1),                            &
  813     &                      ad_zeta(:,:,l1),                            &
  814     &                      nl_zeta(:,:,l1))
  815
  816
  817
  818        fac1=1.0_r8
  819        fac2=zlanczos_inv(ivec,ivec)
  820 
  821        CALL state_addition (ng, tile,                                  &
  822     &                       lbi, ubi, lbj, ubj, lbij, ubij,            &
  823     &                       l1, l1, l1, fac1, fac2,                    &
  824# ifdef MASKING
  825     &                       rmask, umask, vmask,                       &
  826# endif
  827# ifdef ADJUST_BOUNDARY
  828#  ifdef SOLVE3D
  829     &                       tl_t_obc, nl_t_obc,                        &
  830     &                       tl_u_obc, nl_u_obc,                        &
  831     &                       tl_v_obc, nl_v_obc,                        &
  832#  endif
  833     &                       tl_ubar_obc, nl_ubar_obc,                  &
  834     &                       tl_vbar_obc, nl_vbar_obc,                  &
  835     &                       tl_zeta_obc, nl_zeta_obc,                  &
  836# endif
  837# ifdef ADJUST_WSTRESS
  838     &                       tl_ustr, nl_ustr,                          &
  839     &                       tl_vstr, nl_vstr,                          &
  840# endif
  841# ifdef SOLVE3D
  842#  ifdef ADJUST_STFLUX
  843     &                       tl_tflux, nl_tflux,                        &
  844#  endif
  845     &                       tl_t, nl_t,                                &
  846     &                       tl_u, nl_u,                                &
  847     &                       tl_v, nl_v,                                &
  848#  if defined WEAK_CONSTRAINT && defined TIME_CONV
  849     &                       tl_ubar, nl_ubar,                          &
  850     &                       tl_vbar, nl_vbar,                          &
  851#  endif
  852# else
  853     &                       tl_ubar, nl_ubar,                          &
  854     &                       tl_vbar, nl_vbar,                          &
  855# endif
  856     &                       tl_zeta, nl_zeta)
  857
  858        DO jvec=ivec+1,ninner
  859          rec=jvec
  860
  861
  862
  863
  864          CALL state_read (ng, tile, model, hss(ng)%IOtype,             &
  865     &                     lbi, ubi, lbj, ubj, lbij, ubij,              &
  866     &                     l2, rec,                                     &
  867     &                     0, hss(ng)%ncid, ncname,                     &
  868# ifdef MASKING
  869     &                     rmask, umask, vmask,                         &
  870# endif
  871# ifdef ADJUST_BOUNDARY
  872#  ifdef SOLVE3D
  873     &                     ad_t_obc, ad_u_obc, ad_v_obc,                &
  874#  endif
  875     &                     ad_ubar_obc, ad_vbar_obc,                    &
  876     &                     ad_zeta_obc,                                 &
  877# endif
  878# ifdef ADJUST_WSTRESS
  879     &                     ad_ustr, ad_vstr,                            &
  880# endif
  881# ifdef SOLVE3D
  882#  ifdef ADJUST_STFLUX
  883     &                     ad_tflux,                                    &
  884#  endif
  885     &                     ad_t, ad_u, ad_v,                            &
  886# else
  887     &                     ad_ubar, ad_vbar,                            &
  888# endif
  889     &                     ad_zeta)
  890
  891
  892
  893
  894          CALL state_product (ng, tile, model,                          &
  895     &                        lbi, ubi, lbj, ubj, lbij, ubij,           &
  896# ifdef MASKING
  897     &                        rmask, umask, vmask,                      &
  898# endif
  899# ifdef ADJUST_BOUNDARY
  900#  ifdef SOLVE3D
  901     &                        ad_t_obc(:,:,:,:,l1,:),                   &
  902     &                        ad_t_obc(:,:,:,:,l2,:),                   &
  903     &                        nl_t_obc(:,:,:,:,l1,:),                   &
  904     &                        ad_u_obc(:,:,:,:,l1),                     &
  905     &                        ad_u_obc(:,:,:,:,l2),                     &
  906     &                        nl_u_obc(:,:,:,:,l1),                     &
  907     &                        ad_v_obc(:,:,:,:,l1),                     &
  908     &                        ad_v_obc(:,:,:,:,l2),                     &
  909     &                        nl_v_obc(:,:,:,:,l1),                     &
  910#  endif
  911     &                        ad_ubar_obc(:,:,:,l1),                    &
  912     &                        ad_ubar_obc(:,:,:,l2),                    &
  913     &                        nl_ubar_obc(:,:,:,l1),                    &
  914     &                        ad_vbar_obc(:,:,:,l1),                    &
  915     &                        ad_vbar_obc(:,:,:,l2),                    &
  916     &                        nl_vbar_obc(:,:,:,l1),                    &
  917     &                        ad_zeta_obc(:,:,:,l1),                    &
  918     &                        ad_zeta_obc(:,:,:,l2),                    &
  919     &                        nl_zeta_obc(:,:,:,l1),                    &
  920# endif
  921# ifdef ADJUST_WSTRESS
  922     &                        ad_ustr(:,:,:,l1),                        &
  923     &                        ad_ustr(:,:,:,l2),                        &
  924     &                        nl_ustr(:,:,:,l1),                        &
  925     &                        ad_vstr(:,:,:,l1),                        &
  926     &                        ad_vstr(:,:,:,l2),                        &
  927     &                        nl_vstr(:,:,:,l1),                        &
  928# endif
  929# ifdef SOLVE3D
  930#  ifdef ADJUST_STFLUX
  931     &                        ad_tflux(:,:,:,l1,:),                     &
  932     &                        ad_tflux(:,:,:,l2,:),                     &
  933     &                        nl_tflux(:,:,:,l1,:),                     &
  934#  endif
  935     &                        ad_t(:,:,:,l1,:),                         &
  936     &                        ad_t(:,:,:,l2,:),                         &
  937     &                        nl_t(:,:,:,l1,:),                         &
  938     &                        ad_u(:,:,:,l1),                           &
  939     &                        ad_u(:,:,:,l2),                           &
  940     &                        nl_u(:,:,:,l1),                           &
  941     &                        ad_v(:,:,:,l1),                           &
  942     &                        ad_v(:,:,:,l2),                           &
  943     &                        nl_v(:,:,:,l1),                           &
  944# else
  945     &                        ad_ubar(:,:,l1),                          &
  946     &                        ad_ubar(:,:,l2),                          &
  947     &                        nl_ubar(:,:,l1),                          &
  948     &                        ad_vbar(:,:,l1),                          &
  949     &                        ad_vbar(:,:,l2),                          &
  950     &                        nl_vbar(:,:,l1),                          &
  951# endif
  952     &                        ad_zeta(:,:,l1),                          &
  953     &                        ad_zeta(:,:,l2),                          &
  954     &                        nl_zeta(:,:,l1))
  955
  956
  957
  958          fac1=1.0_r8
  959          fac2=2.0_r8*zlanczos_inv(ivec,jvec)
  960 
  961          CALL state_addition (ng, tile,                                &
  962     &                         lbi, ubi, lbj, ubj, lbij, ubij,          &
  963     &                         l1, l1, l1, fac1, fac2,                  &
  964# ifdef MASKING
  965     &                         rmask, umask, vmask,                     &
  966# endif
  967# ifdef ADJUST_BOUNDARY
  968#  ifdef SOLVE3D
  969     &                         tl_t_obc, nl_t_obc,                      &
  970     &                         tl_u_obc, nl_u_obc,                      &
  971     &                         tl_v_obc, nl_v_obc,                      &
  972#  endif
  973     &                         tl_ubar_obc, nl_ubar_obc,                &
  974     &                         tl_vbar_obc, nl_vbar_obc,                &
  975     &                         tl_zeta_obc, nl_zeta_obc,                &
  976# endif
  977# ifdef ADJUST_WSTRESS
  978     &                         tl_ustr, nl_ustr,                        &
  979     &                         tl_vstr, nl_vstr,                        &
  980# endif
  981# ifdef SOLVE3D
  982#  ifdef ADJUST_STFLUX
  983     &                         tl_tflux, nl_tflux,                      &
  984#  endif
  985     &                         tl_t, nl_t,                              &
  986     &                         tl_u, nl_u,                              &
  987     &                         tl_v, nl_v,                              &
  988#  if defined WEAK_CONSTRAINT && defined TIME_CONV
  989     &                         tl_ubar, nl_ubar,                        &
  990     &                         tl_vbar, nl_vbar,                        &
  991#  endif
  992# else
  993     &                         tl_ubar, nl_ubar,                        &
  994     &                         tl_vbar, nl_vbar,                        &
  995# endif
  996     &                         tl_zeta, nl_zeta)
  997        END DO
  998      END DO
  999
 1000
 1001
 1002
 1003      CALL analysis_var (ng, tile, model,                               &
 1004     &                   lbi, ubi, lbj, ubj, lbij, ubij,                &
 1005     &                   l1, lweak,                                     &
 1006# ifdef MASKING
 1007     &                   rmask, umask, vmask,                           &
 1008# endif
 1009# ifdef ADJUST_BOUNDARY
 1010#  ifdef SOLVE3D
 1011     &                   t_obc_std, u_obc_std, v_obc_std,               &
 1012#  endif
 1013     &                   ubar_obc_std, vbar_obc_std, zeta_obc_std,      &
 1014# endif
 1015# ifdef ADJUST_WSTRESS
 1016     &                   sustr_std, svstr_std,                          &
 1017# endif
 1018# ifdef SOLVE3D
 1019#  ifdef ADJUST_STFLUX
 1020     &                   stflx_std,                                     &
 1021#  endif
 1022     &                   t_std, u_std, v_std,                           &
 1023# else
 1024     &                   ubar_std, vbar_std,                            &
 1025# endif
 1026     &                   zeta_std,                                      &
 1027# ifdef ADJUST_BOUNDARY
 1028#  ifdef SOLVE3D
 1029     &                   tl_t_obc, tl_u_obc, tl_v_obc,                  &
 1030#  endif
 1031     &                   tl_ubar_obc, tl_vbar_obc, tl_zeta_obc,         &
 1032# endif
 1033# ifdef ADJUST_WSTRESS
 1034     &                   tl_ustr, tl_vstr,                              &
 1035# endif
 1036# ifdef SOLVE3D
 1037#  ifdef ADJUST_STFLUX
 1038     &                   tl_tflux,                                      &
 1039#  endif
 1040     &                   tl_t, tl_u, tl_v,                              &
 1041# else
 1042     &                   tl_ubar, tl_vbar,                              &
 1043# endif
 1044     &                   tl_zeta)
 1045
 1046      RETURN