129
  130
  133
  134
  135
  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
  140 
  141#ifdef ASSUMED_SHAPE
  142# ifdef MASKING
  143      real(r8), intent(in) :: umask(LBi:,LBj:)
  144      real(r8), intent(in) :: vmask(LBi:,LBj:)
  145# endif
  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:,:)
  152 
  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:)
  160# endif
  161# ifdef ATM_PRESS
  162      real(r8), intent(in) :: Pair(LBi:,LBj:)
  163# endif
  164# ifdef DIAGNOSTICS_UV
  165
  166
  167# endif
  168      real(r8), intent(inout) :: tl_ru(LBi:,LBj:,0:,:)
  169      real(r8), intent(inout) :: tl_rv(LBi:,LBj:,0:,:)
  170#else
  171# ifdef MASKING
  172      real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
  173      real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
  174# endif
  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))
  181 
  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)
  189# endif
  190# ifdef ATM_PRESS
  191      real(r8), intent(in) :: Pair(LBi:UBi,LBj:UBj)
  192# endif
  193# ifdef DIAGNOSTICS_UV
  194
  195
  196# endif
  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)
  199#endif
  200
  201
  202
  203      integer :: i, j, k
  204 
  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
  208 
  209      real(r8) :: GRho, GRho0, HalfGRho
  210      real(r8) :: cff, cff1, cff2
  211      real(r8) :: tl_cff, tl_cff1, tl_cff2
  212#ifdef ATM_PRESS
  213      real(r8) :: OneAtm, fac
  214#endif
  215      real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: P
  216 
  217      real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: tl_P
  218 
  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
  223 
  224      real(r8), dimension(IminS:ImaxS,0:N(ng)) :: tl_dR
  225      real(r8), dimension(IminS:ImaxS,0:N(ng)) :: tl_dZ
  226 
  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
  231 
  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
  236 
  237#include "set_bounds.h"
  238
  239
  240
  241
  242
  244      grho0=1000.0_r8*grho
  245      halfgrho=0.5_r8*grho
  246#ifdef ATM_PRESS
  247      oneatm=1013.25_r8                  
  249#endif
  250
  251      DO j=jstrv-1,jend
  253          DO i=istru-1,iend
  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)
  258          END DO
  259        END DO
  260        DO i=istru-1,iend
  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)
 
  265          dr(i,0)=dr(i,1)
  266          tl_dr(i,0)=tl_dr(i,1)
  267          dz(i,0)=dz(i,1)
  268          tl_dz(i,0)=tl_dz(i,1)
  269        END DO
  271          DO i=istru-1,iend
  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))
  275            dr1(i,k)=dr(i,k)
  276            IF (cff.gt.eps) THEN
  277              dr(i,k)=cff/(dr(i,k)+dr(i,k-1))
  278              tl_dr(i,k)=(tl_cff-dr(i,k)*(tl_dr(i,k)+tl_dr(i,k-1)))/    &
  279     &                   (dr1(i,k)+dr(i,k-1))
  280            ELSE
  281              dr(i,k)=0.0_r8
  282              tl_dr(i,k)=0.0_r8
  283            END IF
  284            dz1(i,k)=dz(i,k)
  285            dz(i,k)=2.0_r8*dz(i,k)*dz(i,k-1)/(dz(i,k)+dz(i,k-1))
  286            tl_dz(i,k)=(2.0_r8*(tl_dz(i,k)*dz(i,k-1)+                   &
  287     &                          dz1(i,k)*tl_dz(i,k-1))-                 &
  288     &                  dz(i,k)*(tl_dz(i,k)+tl_dz(i,k-1)))/             &
  289     &                 (dz1(i,k)+dz(i,k-1))
  290          END DO
  291        END DO
  292        DO i=istru-1,iend
  293          cff1=1.0_r8/(z_r(i,j,
n(ng))-z_r(i,j,
n(ng)-1))
 
  294          tl_cff1=-cff1*cff1*(tl_z_r(i,j,
n(ng))-tl_z_r(i,j,
n(ng)-1))
 
  295          cff2=0.5_r8*(rho(i,j,
n(ng))-rho(i,j,
n(ng)-1))*                &
 
  296     &         (z_w(i,j,
n(ng))-z_r(i,j,
n(ng)))*cff1
 
  297          tl_cff2=0.5_r8*((tl_rho(i,j,
n(ng))-tl_rho(i,j,
n(ng)-1))*      &
 
  298     &                    (z_w(i,j,
n(ng))-z_r(i,j,
n(ng)))*cff1+         &
 
  299     &                    (rho(i,j,
n(ng))-rho(i,j,
n(ng)-1))*            &
 
  300     &                    ((tl_z_w(i,j,
n(ng))-tl_z_r(i,j,
n(ng)))*cff1+  &
 
  301     &                     (z_w(i,j,
n(ng))-z_r(i,j,
n(ng)))*tl_cff1))
 
  302          p(i,j,
n(ng))=
g*z_w(i,j,
n(ng))+                                &
 
  303#ifdef ATM_PRESS
  304     &                 fac*(pair(i,j)-oneatm)+                          &
  305#endif
  306     &                 grho*(rho(i,j,
n(ng))+cff2)*                      &
 
  307     &                 (z_w(i,j,
n(ng))-z_r(i,j,
n(ng)))
 
  308          tl_p(i,j,
n(ng))=
g*tl_z_w(i,j,
n(ng))+                          &
 
  309     &                    grho*((tl_rho(i,j,
n(ng))+tl_cff2)*            &
 
  310     &                          (z_w(i,j,
n(ng))-z_r(i,j,
n(ng)))+        &
 
  311     &                          (rho(i,j,
n(ng))+cff2)*                  &
 
  312     &                          (tl_z_w(i,j,
n(ng))-tl_z_r(i,j,
n(ng))))
 
  313#ifdef TIDE_GENERATING_FORCES
  314          p(i,j,
n(ng))=p(i,j,
n(ng))-
g*eq_tide(i,j)
 
  315          tl_p(i,j,
n(ng))=tl_p(i,j,
n(ng))-
g*tl_eq_tide(i,j)
 
  316#endif
  317        END DO
  319          DO i=istru-1,iend
  320            cff=halfgrho*((rho(i,j,k+1)+rho(i,j,k))*                    &
  321     &                    (z_r(i,j,k+1)-z_r(i,j,k))-                    &
  322     &                    onefifth*                                     &
  323     &                    ((dr(i,k+1)-dr(i,k))*                         &
  324     &                     (z_r(i,j,k+1)-z_r(i,j,k)-                    &
  325     &                      onetwelfth*                                 &
  326     &                      (dz(i,k+1)+dz(i,k)))-                       &
  327     &                     (dz(i,k+1)-dz(i,k))*                         &
  328     &                     (rho(i,j,k+1)-rho(i,j,k)-                    &
  329     &                      onetwelfth*                                 &
  330     &                      (dr(i,k+1)+dr(i,k)))))
  331            tl_cff=halfgrho*((tl_rho(i,j,k+1)+tl_rho(i,j,k))*           &
  332     &                       (z_r(i,j,k+1)-z_r(i,j,k))+                 &
  333     &                       (rho(i,j,k+1)+rho(i,j,k))*                 &
  334     &                       (tl_z_r(i,j,k+1)-tl_z_r(i,j,k))-           &
  335     &                       onefifth*                                  &
  336     &                       ((tl_dr(i,k+1)-tl_dr(i,k))*                &
  337     &                        (z_r(i,j,k+1)-z_r(i,j,k)-                 &
  338     &                         onetwelfth*                              &
  339     &                         (dz(i,k+1)+dz(i,k)))+                    &
  340     &                        (dr(i,k+1)-dr(i,k))*                      &
  341     &                        (tl_z_r(i,j,k+1)-tl_z_r(i,j,k)-           &
  342     &                         onetwelfth*                              &
  343     &                         (tl_dz(i,k+1)+tl_dz(i,k)))-              &
  344     &                        (tl_dz(i,k+1)-tl_dz(i,k))*                &
  345     &                        (rho(i,j,k+1)-rho(i,j,k)-                 &
  346     &                         onetwelfth*                              &
  347     &                         (dr(i,k+1)+dr(i,k)))-                    &
  348     &                        (dz(i,k+1)-dz(i,k))*                      &
  349     &                        (tl_rho(i,j,k+1)-tl_rho(i,j,k)-           &
  350     &                         onetwelfth*                              &
  351     &                         (tl_dr(i,k+1)+tl_dr(i,k)))))
  352            p(i,j,k)=p(i,j,k+1)+cff
  353            tl_p(i,j,k)=tl_p(i,j,k+1)+tl_cff
  354          END DO
  355        END DO
  356      END DO
  357
  358
  359
  360
  361
  363        DO j=jstr,jend
  364          DO i=istru-1,iend+1
  365            aux(i,j)=z_r(i,j,k)-z_r(i-1,j,k)
  366            tl_aux(i,j)=tl_z_r(i,j,k)-tl_z_r(i-1,j,k)
  367#ifdef MASKING
  368            aux(i,j)=aux(i,j)*umask(i,j)
  369            tl_aux(i,j)=tl_aux(i,j)*umask(i,j)
  370#endif
  371            fc(i,j)=rho(i,j,k)-rho(i-1,j,k)
  372            tl_fc(i,j)=tl_rho(i,j,k)-tl_rho(i-1,j,k)
  373#ifdef MASKING
  374            fc(i,j)=fc(i,j)*umask(i,j)
  375            tl_fc(i,j)=tl_fc(i,j)*umask(i,j)
  376#endif
  377          END DO
  378        END DO
  379
  380        DO j=jstr,jend
  381          DO i=istru-1,iend
  382            cff=2.0_r8*aux(i,j)*aux(i+1,j)
  383            tl_cff=2.0_r8*(tl_aux(i,j)*aux(i+1,j)+                      &
  384     &                     aux(i,j)*tl_aux(i+1,j))
  385            IF (cff.gt.eps) THEN
  386              cff1=1.0_r8/(aux(i,j)+aux(i+1,j))
  387              tl_cff1=-cff1*cff1*(tl_aux(i,j)+tl_aux(i+1,j))
  388              dzx(i,j)=cff*cff1
  389              tl_dzx(i,j)=tl_cff*cff1+cff*tl_cff1
  390            ELSE
  391              dzx(i,j)=0.0_r8
  392              tl_dzx(i,j)=0.0_r8
  393            END IF
  394            cff1=2.0_r8*fc(i,j)*fc(i+1,j)
  395            tl_cff1=2.0_r8*(tl_fc(i,j)*fc(i+1,j)+                       &
  396     &                      fc(i,j)*tl_fc(i+1,j))
  397            IF (cff1.gt.eps) THEN
  398              cff2=1.0_r8/(fc(i,j)+fc(i+1,j))
  399              tl_cff2=-cff2*cff2*(tl_fc(i,j)+tl_fc(i+1,j))
  400              drx(i,j)=cff1*cff2
  401              tl_drx(i,j)=tl_cff1*cff2+cff1*tl_cff2
  402            ELSE
  403              drx(i,j)=0.0_r8
  404              tl_drx(i,j)=0.0_r8
  405            END IF
  406          END DO
  407        END DO
  408
  409        DO j=jstr,jend
  410          DO i=istru,iend
  411
  412
  413
  414
  415
  416
  417
  418
  419
  420
  421
  422
  423
  424
  425
  426
  427            tl_ru(i,j,k,nrhs)=on_u(i,j)*0.5_r8*                         &
  428     &                        ((tl_hz(i,j,k)+tl_hz(i-1,j,k))*           &
  429     &                         (p(i-1,j,k)-p(i,j,k)-                    &
  430     &                          halfgrho*                               &
  431     &                          ((rho(i,j,k)+rho(i-1,j,k))*             &
  432     &                           (z_r(i,j,k)-z_r(i-1,j,k))-             &
  433     &                            onefifth*                             &
  434     &                            ((drx(i,j)-drx(i-1,j))*               &
  435     &                             (z_r(i,j,k)-z_r(i-1,j,k)-            &
  436     &                              onetwelfth*                         &
  437     &                              (dzx(i,j)+dzx(i-1,j)))-             &
  438     &                             (dzx(i,j)-dzx(i-1,j))*               &
  439     &                             (rho(i,j,k)-rho(i-1,j,k)-            &
  440     &                              onetwelfth*                         &
  441     &                              (drx(i,j)+drx(i-1,j))))))+          &
  442     &                         (hz(i,j,k)+hz(i-1,j,k))*                 &
  443     &                         (tl_p(i-1,j,k)-tl_p(i,j,k)-              &
  444     &                          halfgrho*                               &
  445     &                          ((tl_rho(i,j,k)+tl_rho(i-1,j,k))*       &
  446     &                           (z_r(i,j,k)-z_r(i-1,j,k))+             &
  447     &                           (rho(i,j,k)+rho(i-1,j,k))*             &
  448     &                           (tl_z_r(i,j,k)-tl_z_r(i-1,j,k))-       &
  449     &                            onefifth*                             &
  450     &                            ((tl_drx(i,j)-tl_drx(i-1,j))*         &
  451     &                             (z_r(i,j,k)-z_r(i-1,j,k)-            &
  452     &                              onetwelfth*                         &
  453     &                              (dzx(i,j)+dzx(i-1,j)))+             &
  454     &                             (drx(i,j)-drx(i-1,j))*               &
  455     &                             (tl_z_r(i,j,k)-tl_z_r(i-1,j,k)-      &
  456     &                              onetwelfth*                         &
  457     &                              (tl_dzx(i,j)+tl_dzx(i-1,j)))-       &
  458     &                             (tl_dzx(i,j)-tl_dzx(i-1,j))*         &
  459     &                             (rho(i,j,k)-rho(i-1,j,k)-            &
  460     &                              onetwelfth*                         &
  461     &                              (drx(i,j)+drx(i-1,j)))-             &
  462     &                             (dzx(i,j)-dzx(i-1,j))*               &
  463     &                             (tl_rho(i,j,k)-tl_rho(i-1,j,k)-      &
  464     &                              onetwelfth*                         &
  465     &                              (tl_drx(i,j)+tl_drx(i-1,j)))))))
  466#ifdef DIAGNOSTICS_UV
  467
  468#endif
  469          END DO
  470        END DO
  471      END DO
  472
  473
  474
  475
  476
  478        DO j=jstrv-1,jend+1
  479          DO i=istr,iend
  480            aux(i,j)=z_r(i,j,k)-z_r(i,j-1,k)
  481            tl_aux(i,j)=tl_z_r(i,j,k)-tl_z_r(i,j-1,k)
  482#ifdef MASKING
  483            aux(i,j)=aux(i,j)*vmask(i,j)
  484            tl_aux(i,j)=tl_aux(i,j)*vmask(i,j)
  485#endif
  486            fc(i,j)=rho(i,j,k)-rho(i,j-1,k)
  487            tl_fc(i,j)=tl_rho(i,j,k)-tl_rho(i,j-1,k)
  488#ifdef MASKING
  489            fc(i,j)=fc(i,j)*vmask(i,j)
  490            tl_fc(i,j)=tl_fc(i,j)*vmask(i,j)
  491#endif
  492          END DO
  493        END DO
  494
  495        DO j=jstrv-1,jend
  496          DO i=istr,iend
  497            cff=2.0_r8*aux(i,j)*aux(i,j+1)
  498            tl_cff=2.0_r8*(tl_aux(i,j)*aux(i,j+1)+                      &
  499     &                     aux(i,j)*tl_aux(i,j+1))
  500            IF (cff.gt.eps) THEN
  501              cff1=1.0_r8/(aux(i,j)+aux(i,j+1))
  502              tl_cff1=-cff1*cff1*(tl_aux(i,j)+tl_aux(i,j+1))
  503              dzx(i,j)=cff*cff1
  504              tl_dzx(i,j)=tl_cff*cff1+cff*tl_cff1
  505            ELSE
  506              dzx(i,j)=0.0_r8
  507              tl_dzx(i,j)=0.0_r8
  508            END IF
  509            cff1=2.0_r8*fc(i,j)*fc(i,j+1)
  510            tl_cff1=2.0_r8*(tl_fc(i,j)*fc(i,j+1)+                       &
  511     &                      fc(i,j)*tl_fc(i,j+1))
  512            IF (cff1.gt.eps) THEN
  513              cff2=1.0_r8/(fc(i,j)+fc(i,j+1))
  514              tl_cff2=-cff2*cff2*(tl_fc(i,j)+tl_fc(i,j+1))
  515              drx(i,j)=cff1*cff2
  516              tl_drx(i,j)=tl_cff1*cff2+cff1*tl_cff2
  517            ELSE
  518              drx(i,j)=0.0_r8
  519              tl_drx(i,j)=0.0_r8
  520            END IF
  521          END DO
  522        END DO
  523
  524        DO j=jstrv,jend
  525          DO i=istr,iend
  526
  527
  528
  529
  530
  531
  532
  533
  534
  535
  536
  537
  538
  539
  540
  541
  542            tl_rv(i,j,k,nrhs)=om_v(i,j)*0.5_r8*                         &
  543     &                        ((tl_hz(i,j,k)+tl_hz(i,j-1,k))*           &
  544     &                         (p(i,j-1,k)-p(i,j,k)-                    &
  545     &                          halfgrho*                               &
  546     &                          ((rho(i,j,k)+rho(i,j-1,k))*             &
  547     &                           (z_r(i,j,k)-z_r(i,j-1,k))-             &
  548     &                           onefifth*                              &
  549     &                           ((drx(i,j)-drx(i,j-1))*                &
  550     &                            (z_r(i,j,k)-z_r(i,j-1,k)-             &
  551     &                             onetwelfth*                          &
  552     &                             (dzx(i,j)+dzx(i,j-1)))-              &
  553     &                            (dzx(i,j)-dzx(i,j-1))*                &
  554     &                            (rho(i,j,k)-rho(i,j-1,k)-             &
  555     &                             onetwelfth*                          &
  556     &                             (drx(i,j)+drx(i,j-1))))))+           &
  557     &                         (hz(i,j,k)+hz(i,j-1,k))*                 &
  558     &                         (tl_p(i,j-1,k)-tl_p(i,j,k)-              &
  559     &                          halfgrho*                               &
  560     &                          ((tl_rho(i,j,k)+tl_rho(i,j-1,k))*       &
  561     &                           (z_r(i,j,k)-z_r(i,j-1,k))+             &
  562     &                           (rho(i,j,k)+rho(i,j-1,k))*             &
  563     &                           (tl_z_r(i,j,k)-tl_z_r(i,j-1,k))-       &
  564     &                           onefifth*                              &
  565     &                           ((tl_drx(i,j)-tl_drx(i,j-1))*          &
  566     &                            (z_r(i,j,k)-z_r(i,j-1,k)-             &
  567     &                             onetwelfth*                          &
  568     &                             (dzx(i,j)+dzx(i,j-1)))+              &
  569     &                            (drx(i,j)-drx(i,j-1))*                &
  570     &                            (tl_z_r(i,j,k)-tl_z_r(i,j-1,k)-       &
  571     &                             onetwelfth*                          &
  572     &                             (tl_dzx(i,j)+tl_dzx(i,j-1)))-        &
  573     &                            (tl_dzx(i,j)-tl_dzx(i,j-1))*          &
  574     &                            (rho(i,j,k)-rho(i,j-1,k)-             &
  575     &                             onetwelfth*                          &
  576     &                             (drx(i,j)+drx(i,j-1)))-              &
  577     &                            (dzx(i,j)-dzx(i,j-1))*                &
  578     &                            (tl_rho(i,j,k)-tl_rho(i,j-1,k)-       &
  579     &                             onetwelfth*                          &
  580     &                             (tl_drx(i,j)+tl_drx(i,j-1)))))))
  581#ifdef DIAGNOSTICS_UV
  582
  583#endif
  584          END DO
  585        END DO
  586      END DO
  587
  588      RETURN