108     &                             LBi, UBi, LBj, UBj,                  &
 
  109     &                             IminS, ImaxS, JminS, JmaxS,          &
 
  119#ifdef TIDE_GENERATING_FORCES
 
  120     &                             eq_tide, tl_eq_tide,                 &
 
  136      integer, 
intent(in) :: ng, tile
 
  137      integer, 
intent(in) :: LBi, UBi, LBj, UBj
 
  138      integer, 
intent(in) :: IminS, ImaxS, JminS, JmaxS
 
  139      integer, 
intent(in) :: nrhs
 
  143      real(r8), 
intent(in) :: umask(LBi:,LBj:)
 
  144      real(r8), 
intent(in) :: vmask(LBi:,LBj:)
 
  146      real(r8), 
intent(in) :: om_v(LBi:,LBj:)
 
  147      real(r8), 
intent(in) :: on_u(LBi:,LBj:)
 
  148      real(r8), 
intent(in) :: Hz(LBi:,LBj:,:)
 
  149      real(r8), 
intent(in) :: z_r(LBi:,LBj:,:)
 
  150      real(r8), 
intent(in) :: z_w(LBi:,LBj:,0:)
 
  151      real(r8), 
intent(in) :: rho(LBi:,LBj:,:)
 
  153      real(r8), 
intent(in) :: tl_Hz(LBi:,LBj:,:)
 
  154      real(r8), 
intent(in) :: tl_z_r(LBi:,LBj:,:)
 
  155      real(r8), 
intent(in) :: tl_z_w(LBi:,LBj:,0:)
 
  156      real(r8), 
intent(in) :: tl_rho(LBi:,LBj:,:)
 
  157# ifdef TIDE_GENERATING_FORCES 
  158      real(r8), 
intent(in) :: eq_tide(LBi:,LBj:)
 
  159      real(r8), 
intent(in) :: tl_eq_tide(LBi:,LBj:)
 
  162      real(r8), 
intent(in) :: Pair(LBi:,LBj:)
 
  164# ifdef DIAGNOSTICS_UV 
  168      real(r8), 
intent(inout) :: tl_ru(LBi:,LBj:,0:,:)
 
  169      real(r8), 
intent(inout) :: tl_rv(LBi:,LBj:,0:,:)
 
  172      real(r8), 
intent(in) :: umask(LBi:UBi,LBj:UBj)
 
  173      real(r8), 
intent(in) :: vmask(LBi:UBi,LBj:UBj)
 
  175      real(r8), 
intent(in) :: om_v(LBi:UBi,LBj:UBj)
 
  176      real(r8), 
intent(in) :: on_u(LBi:UBi,LBj:UBj)
 
  177      real(r8), 
intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
 
  178      real(r8), 
intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
 
  179      real(r8), 
intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
 
  180      real(r8), 
intent(in) :: rho(LBi:UBi,LBj:UBj,N(ng))
 
  182      real(r8), 
intent(in) :: tl_Hz(LBi:UBi,LBj:UBj,N(ng))
 
  183      real(r8), 
intent(in) :: tl_z_r(LBi:UBi,LBj:UBj,N(ng))
 
  184      real(r8), 
intent(in) :: tl_z_w(LBi:UBi,LBj:UBj,0:N(ng))
 
  185      real(r8), 
intent(in) :: tl_rho(LBi:UBi,LBj:UBj,N(ng))
 
  186# ifdef TIDE_GENERATING_FORCES 
  187      real(r8), 
intent(in) :: eq_tide(LBi:UBi,LBj:UBj)
 
  188      real(r8), 
intent(in) :: tl_eq_tide(LBi:UBi,LBj:UBj)
 
  191      real(r8), 
intent(in) :: Pair(LBi:UBi,LBj:UBj)
 
  193# ifdef DIAGNOSTICS_UV 
  197      real(r8), 
intent(inout) :: tl_ru(LBi:UBi,LBj:UBj,0:N(ng),2)
 
  198      real(r8), 
intent(inout) :: tl_rv(LBi:UBi,LBj:UBj,0:N(ng),2)
 
  205      real(r8), 
parameter :: OneFifth = 0.2_r8
 
  206      real(r8), 
parameter :: OneTwelfth = 1.0_r8/12.0_r8
 
  207      real(r8), 
parameter :: eps = 1.0e-10_r8
 
  209      real(r8) :: GRho, GRho0, HalfGRho
 
  210      real(r8) :: cff, cff1, cff2
 
  211      real(r8) :: tl_cff, tl_cff1, tl_cff2
 
  213      real(r8) :: OneAtm, fac
 
  215      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: P
 
  217      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: tl_P
 
  219      real(r8), 
dimension(IminS:ImaxS,0:N(ng)) :: dR
 
  220      real(r8), 
dimension(IminS:ImaxS,0:N(ng)) :: dR1
 
  221      real(r8), 
dimension(IminS:ImaxS,0:N(ng)) :: dZ
 
  222      real(r8), 
dimension(IminS:ImaxS,0:N(ng)) :: dZ1
 
  224      real(r8), 
dimension(IminS:ImaxS,0:N(ng)) :: tl_dR
 
  225      real(r8), 
dimension(IminS:ImaxS,0:N(ng)) :: tl_dZ
 
  227      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS) :: FC
 
  228      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS) :: aux
 
  229      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS) :: dRx
 
  230      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS) :: dZx
 
  232      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS) :: tl_FC
 
  233      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS) :: tl_aux
 
  234      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS) :: tl_dRx
 
  235      real(r8), 
dimension(IminS:ImaxS,JminS:JmaxS) :: tl_dZx
 
  237#include "set_bounds.h" 
  254            dr(i,k)=rho(i,j,k+1)-rho(i,j,k)
 
  255            tl_dr(i,k)=tl_rho(i,j,k+1)-tl_rho(i,j,k)
 
  256            dz(i,k)=z_r(i,j,k+1)-z_r(i,j,k)
 
  257            tl_dz(i,k)=tl_z_r(i,j,k+1)-tl_z_r(i,j,k)
 
  261          dr(i,n(ng))=dr(i,n(ng)-1)
 
  262          tl_dr(i,n(ng))=tl_dr(i,n(ng)-1)
 
  263          dz(i,n(ng))=dz(i,n(ng)-1)
 
  264          tl_dz(i,n(ng))=tl_dz(i,n(ng)-1)
 
  266          tl_dr(i,0)=tl_dr(i,1)
 
  268          tl_dz(i,0)=tl_dz(i,1)
 
  272            cff=2.0_r8*dr(i,k)*dr(i,k-1)
 
  273            tl_cff=2.0_r8*(tl_dr(i,k)*dr(i,k-1)+                        &
 
  274     &                     dr(i,k)*tl_dr(i,k-1))-                       &
 
  280              dr(i,k)=cff/(dr(i,k)+dr(i,k-1))
 
  281              tl_dr(i,k)=(tl_cff-dr(i,k)*(tl_dr(i,k)+tl_dr(i,k-1)))/    &
 
  282     &                   (dr1(i,k)+dr(i,k-1))+                          &
 
  291            dz(i,k)=2.0_r8*dz(i,k)*dz(i,k-1)/(dz(i,k)+dz(i,k-1))
 
  292            tl_dz(i,k)=(2.0_r8*(tl_dz(i,k)*dz(i,k-1)+                   &
 
  293     &                          dz1(i,k)*tl_dz(i,k-1))-                 &
 
  294     &                  dz(i,k)*(tl_dz(i,k)+tl_dz(i,k-1)))/             &
 
  295     &                 (dz1(i,k)+dz(i,k-1))
 
  299          cff1=1.0_r8/(z_r(i,j,n(ng))-z_r(i,j,n(ng)-1))
 
  300          tl_cff1=-cff1*cff1*(tl_z_r(i,j,n(ng))-tl_z_r(i,j,n(ng)-1))+   &
 
  304          cff2=0.5_r8*(rho(i,j,n(ng))-rho(i,j,n(ng)-1))*                &
 
  305     &         (z_w(i,j,n(ng))-z_r(i,j,n(ng)))*cff1
 
  306          tl_cff2=0.5_r8*((tl_rho(i,j,n(ng))-tl_rho(i,j,n(ng)-1))*      &
 
  307     &                    (z_w(i,j,n(ng))-z_r(i,j,n(ng)))*cff1+         &
 
  308     &                    (rho(i,j,n(ng))-rho(i,j,n(ng)-1))*            &
 
  309     &                    ((tl_z_w(i,j,n(ng))-tl_z_r(i,j,n(ng)))*cff1+  &
 
  310     &                     (z_w(i,j,n(ng))-z_r(i,j,n(ng)))*tl_cff1))-   &
 
  314          p(i,j,n(ng))=
g*z_w(i,j,n(ng))+                                &
 
  316     &                 fac*(pair(i,j)-oneatm)+                          &
 
  318     &                 grho*(rho(i,j,n(ng))+cff2)*                      &
 
  319     &                 (z_w(i,j,n(ng))-z_r(i,j,n(ng)))
 
  320          tl_p(i,j,n(ng))=
g*tl_z_w(i,j,n(ng))+                          &
 
  321     &                    grho*((tl_rho(i,j,n(ng))+tl_cff2)*            &
 
  322     &                          (z_w(i,j,n(ng))-z_r(i,j,n(ng)))+        &
 
  323     &                          (rho(i,j,n(ng))+cff2)*                  &
 
  324     &                          (tl_z_w(i,j,n(ng))-tl_z_r(i,j,n(ng))))- &
 
  326     &                    grho*(rho(i,j,n(ng))+cff2)*                   &
 
  327     &                    (z_w(i,j,n(ng))-z_r(i,j,n(ng)))
 
  329#ifdef TIDE_GENERATING_FORCES 
  330          p(i,j,n(ng))=p(i,j,n(ng))-
g*eq_tide(i,j)
 
  331          tl_p(i,j,n(ng))=tl_p(i,j,n(ng))-
g*tl_eq_tide(i,j)
 
  336            cff=halfgrho*((rho(i,j,k+1)+rho(i,j,k))*                    &
 
  337     &                    (z_r(i,j,k+1)-z_r(i,j,k))-                    &
 
  339     &                    ((dr(i,k+1)-dr(i,k))*                         &
 
  340     &                     (z_r(i,j,k+1)-z_r(i,j,k)-                    &
 
  342     &                      (dz(i,k+1)+dz(i,k)))-                       &
 
  343     &                     (dz(i,k+1)-dz(i,k))*                         &
 
  344     &                     (rho(i,j,k+1)-rho(i,j,k)-                    &
 
  346     &                      (dr(i,k+1)+dr(i,k)))))
 
  347            tl_cff=halfgrho*((tl_rho(i,j,k+1)+tl_rho(i,j,k))*           &
 
  348     &                       (z_r(i,j,k+1)-z_r(i,j,k))+                 &
 
  349     &                       (rho(i,j,k+1)+rho(i,j,k))*                 &
 
  350     &                       (tl_z_r(i,j,k+1)-tl_z_r(i,j,k))-           &
 
  352     &                       ((tl_dr(i,k+1)-tl_dr(i,k))*                &
 
  353     &                        (z_r(i,j,k+1)-z_r(i,j,k)-                 &
 
  355     &                         (dz(i,k+1)+dz(i,k)))+                    &
 
  356     &                        (dr(i,k+1)-dr(i,k))*                      &
 
  357     &                        (tl_z_r(i,j,k+1)-tl_z_r(i,j,k)-           &
 
  359     &                         (tl_dz(i,k+1)+tl_dz(i,k)))-              &
 
  360     &                        (tl_dz(i,k+1)-tl_dz(i,k))*                &
 
  361     &                        (rho(i,j,k+1)-rho(i,j,k)-                 &
 
  363     &                         (dr(i,k+1)+dr(i,k)))-                    &
 
  364     &                        (dz(i,k+1)-dz(i,k))*                      &
 
  365     &                        (tl_rho(i,j,k+1)-tl_rho(i,j,k)-           &
 
  367     &                         (tl_dr(i,k+1)+tl_dr(i,k)))))-            &
 
  371            p(i,j,k)=p(i,j,k+1)+cff
 
  372            tl_p(i,j,k)=tl_p(i,j,k+1)+tl_cff
 
  384            aux(i,j)=z_r(i,j,k)-z_r(i-1,j,k)
 
  385            tl_aux(i,j)=tl_z_r(i,j,k)-tl_z_r(i-1,j,k)
 
  387            aux(i,j)=aux(i,j)*umask(i,j)
 
  388            tl_aux(i,j)=tl_aux(i,j)*umask(i,j)
 
  390            fc(i,j)=rho(i,j,k)-rho(i-1,j,k)
 
  391            tl_fc(i,j)=tl_rho(i,j,k)-tl_rho(i-1,j,k)
 
  393            fc(i,j)=fc(i,j)*umask(i,j)
 
  394            tl_fc(i,j)=tl_fc(i,j)*umask(i,j)
 
  401            cff=2.0_r8*aux(i,j)*aux(i+1,j)
 
  402            tl_cff=2.0_r8*(tl_aux(i,j)*aux(i+1,j)+                      &
 
  403     &                     aux(i,j)*tl_aux(i+1,j))-                     &
 
  408              cff1=1.0_r8/(aux(i,j)+aux(i+1,j))
 
  409              tl_cff1=-cff1*cff1*(tl_aux(i,j)+tl_aux(i+1,j))+           &
 
  414              tl_dzx(i,j)=tl_cff*cff1+cff*tl_cff1-                      &
 
  422            cff1=2.0_r8*fc(i,j)*fc(i+1,j)
 
  423            tl_cff1=2.0_r8*(tl_fc(i,j)*fc(i+1,j)+                       &
 
  424     &                      fc(i,j)*tl_fc(i+1,j))-                      &
 
  428            IF (cff1.gt.eps) 
THEN 
  429              cff2=1.0_r8/(fc(i,j)+fc(i+1,j))
 
  430              tl_cff2=-cff2*cff2*(tl_fc(i,j)+tl_fc(i+1,j))+             &
 
  435              tl_drx(i,j)=tl_cff1*cff2+cff1*tl_cff2-                    &
 
  464            tl_ru(i,j,k,nrhs)=on_u(i,j)*0.5_r8*                         &
 
  465     &                        ((tl_hz(i,j,k)+tl_hz(i-1,j,k))*           &
 
  466     &                         (p(i-1,j,k)-p(i,j,k)-                    &
 
  468     &                          ((rho(i,j,k)+rho(i-1,j,k))*             &
 
  469     &                           (z_r(i,j,k)-z_r(i-1,j,k))-             &
 
  471     &                            ((drx(i,j)-drx(i-1,j))*               &
 
  472     &                             (z_r(i,j,k)-z_r(i-1,j,k)-            &
 
  474     &                              (dzx(i,j)+dzx(i-1,j)))-             &
 
  475     &                             (dzx(i,j)-dzx(i-1,j))*               &
 
  476     &                             (rho(i,j,k)-rho(i-1,j,k)-            &
 
  478     &                              (drx(i,j)+drx(i-1,j))))))+          &
 
  479     &                         (hz(i,j,k)+hz(i-1,j,k))*                 &
 
  480     &                         (tl_p(i-1,j,k)-tl_p(i,j,k)-              &
 
  482     &                          ((tl_rho(i,j,k)+tl_rho(i-1,j,k))*       &
 
  483     &                           (z_r(i,j,k)-z_r(i-1,j,k))+             &
 
  484     &                           (rho(i,j,k)+rho(i-1,j,k))*             &
 
  485     &                           (tl_z_r(i,j,k)-tl_z_r(i-1,j,k))-       &
 
  487     &                            ((tl_drx(i,j)-tl_drx(i-1,j))*         &
 
  488     &                             (z_r(i,j,k)-z_r(i-1,j,k)-            &
 
  490     &                              (dzx(i,j)+dzx(i-1,j)))+             &
 
  491     &                             (drx(i,j)-drx(i-1,j))*               &
 
  492     &                             (tl_z_r(i,j,k)-tl_z_r(i-1,j,k)-      &
 
  494     &                              (tl_dzx(i,j)+tl_dzx(i-1,j)))-       &
 
  495     &                             (tl_dzx(i,j)-tl_dzx(i-1,j))*         &
 
  496     &                             (rho(i,j,k)-rho(i-1,j,k)-            &
 
  498     &                              (drx(i,j)+drx(i-1,j)))-             &
 
  499     &                             (dzx(i,j)-dzx(i-1,j))*               &
 
  500     &                             (tl_rho(i,j,k)-tl_rho(i-1,j,k)-      &
 
  502     &                              (tl_drx(i,j)+tl_drx(i-1,j)))))))-   &
 
  504     &                        on_u(i,j)*0.5_r8*                         &
 
  505     &                        (hz(i,j,k)+hz(i-1,j,k))*                  &
 
  506     &                        (p(i-1,j,k)-p(i,j,k)-                     &
 
  508     &                         ((rho(i,j,k)+rho(i-1,j,k))*              &
 
  509     &                          (z_r(i,j,k)-z_r(i-1,j,k))-              &
 
  511     &                           ((drx(i,j)-drx(i-1,j))*                &
 
  512     &                            (z_r(i,j,k)-z_r(i-1,j,k)-             &
 
  514     &                             (dzx(i,j)+dzx(i-1,j)))-              &
 
  515     &                            (dzx(i,j)-dzx(i-1,j))*                &
 
  516     &                            (rho(i,j,k)-rho(i-1,j,k)-             &
 
  518     &                             (drx(i,j)+drx(i-1,j))))))
 
  534            aux(i,j)=z_r(i,j,k)-z_r(i,j-1,k)
 
  535            tl_aux(i,j)=tl_z_r(i,j,k)-tl_z_r(i,j-1,k)
 
  537            aux(i,j)=aux(i,j)*vmask(i,j)
 
  538            tl_aux(i,j)=tl_aux(i,j)*vmask(i,j)
 
  540            fc(i,j)=rho(i,j,k)-rho(i,j-1,k)
 
  541            tl_fc(i,j)=tl_rho(i,j,k)-tl_rho(i,j-1,k)
 
  543            fc(i,j)=fc(i,j)*vmask(i,j)
 
  544            tl_fc(i,j)=tl_fc(i,j)*vmask(i,j)
 
  551            cff=2.0_r8*aux(i,j)*aux(i,j+1)
 
  552            tl_cff=2.0_r8*(tl_aux(i,j)*aux(i,j+1)+                      &
 
  553     &                     aux(i,j)*tl_aux(i,j+1))-                     &
 
  558              cff1=1.0_r8/(aux(i,j)+aux(i,j+1))
 
  559              tl_cff1=-cff1*cff1*(tl_aux(i,j)+tl_aux(i,j+1))+           &
 
  564              tl_dzx(i,j)=tl_cff*cff1+cff*tl_cff1-                      &
 
  572            cff1=2.0_r8*fc(i,j)*fc(i,j+1)
 
  573            tl_cff1=2.0_r8*(tl_fc(i,j)*fc(i,j+1)+                       &
 
  574     &                      fc(i,j)*tl_fc(i,j+1))-                      &
 
  578            IF (cff1.gt.eps) 
THEN 
  579              cff2=1.0_r8/(fc(i,j)+fc(i,j+1))
 
  580              tl_cff2=-cff2*cff2*(tl_fc(i,j)+tl_fc(i,j+1))+             &
 
  585              tl_drx(i,j)=tl_cff1*cff2+cff1*tl_cff2-                    &
 
  614            tl_rv(i,j,k,nrhs)=om_v(i,j)*0.5_r8*                         &
 
  615     &                        ((tl_hz(i,j,k)+tl_hz(i,j-1,k))*           &
 
  616     &                         (p(i,j-1,k)-p(i,j,k)-                    &
 
  618     &                          ((rho(i,j,k)+rho(i,j-1,k))*             &
 
  619     &                           (z_r(i,j,k)-z_r(i,j-1,k))-             &
 
  621     &                           ((drx(i,j)-drx(i,j-1))*                &
 
  622     &                            (z_r(i,j,k)-z_r(i,j-1,k)-             &
 
  624     &                             (dzx(i,j)+dzx(i,j-1)))-              &
 
  625     &                            (dzx(i,j)-dzx(i,j-1))*                &
 
  626     &                            (rho(i,j,k)-rho(i,j-1,k)-             &
 
  628     &                             (drx(i,j)+drx(i,j-1))))))+           &
 
  629     &                         (hz(i,j,k)+hz(i,j-1,k))*                 &
 
  630     &                         (tl_p(i,j-1,k)-tl_p(i,j,k)-              &
 
  632     &                          ((tl_rho(i,j,k)+tl_rho(i,j-1,k))*       &
 
  633     &                           (z_r(i,j,k)-z_r(i,j-1,k))+             &
 
  634     &                           (rho(i,j,k)+rho(i,j-1,k))*             &
 
  635     &                           (tl_z_r(i,j,k)-tl_z_r(i,j-1,k))-       &
 
  637     &                           ((tl_drx(i,j)-tl_drx(i,j-1))*          &
 
  638     &                            (z_r(i,j,k)-z_r(i,j-1,k)-             &
 
  640     &                             (dzx(i,j)+dzx(i,j-1)))+              &
 
  641     &                            (drx(i,j)-drx(i,j-1))*                &
 
  642     &                            (tl_z_r(i,j,k)-tl_z_r(i,j-1,k)-       &
 
  644     &                             (tl_dzx(i,j)+tl_dzx(i,j-1)))-        &
 
  645     &                            (tl_dzx(i,j)-tl_dzx(i,j-1))*          &
 
  646     &                            (rho(i,j,k)-rho(i,j-1,k)-             &
 
  648     &                             (drx(i,j)+drx(i,j-1)))-              &
 
  649     &                            (dzx(i,j)-dzx(i,j-1))*                &
 
  650     &                            (tl_rho(i,j,k)-tl_rho(i,j-1,k)-       &
 
  652     &                             (tl_drx(i,j)+tl_drx(i,j-1)))))))-    &
 
  654     &                        om_v(i,j)*0.5_r8*                         &
 
  655     &                        (hz(i,j,k)+hz(i,j-1,k))*                  &
 
  656     &                        (p(i,j-1,k)-p(i,j,k)-                     &
 
  658     &                         ((rho(i,j,k)+rho(i,j-1,k))*              &
 
  659     &                          (z_r(i,j,k)-z_r(i,j-1,k))-              &
 
  661     &                           ((drx(i,j)-drx(i,j-1))*                &
 
  662     &                            (z_r(i,j,k)-z_r(i,j-1,k)-             &
 
  664     &                             (dzx(i,j)+dzx(i,j-1)))-              &
 
  665     &                            (dzx(i,j)-dzx(i,j-1))*                &
 
  666     &                            (rho(i,j,k)-rho(i,j-1,k)-             &
 
  668     &                             (drx(i,j)+drx(i,j-1))))))