115 & LBi, UBi, LBj, UBj, &
116 & IminS, ImaxS, JminS, JmaxS, &
122 & om_p, om_r, on_p, on_r, &
123 & pm, pmon_p, pmon_r, &
124 & pn, pnom_p, pnom_r, &
125 & visc4_p, visc4_r, &
126!!#ifdef DIAGNOSTICS_UV
127!! & DiaRUfrc, DiaRVfrc, &
128!! & DiaU3wrk, DiaV3wrk, &
131 & ad_rufrc, ad_rvfrc, ad_u, ad_v)
139 integer,
intent(in) :: ng, tile
140 integer,
intent(in) :: LBi, UBi, LBj, UBj
141 integer,
intent(in) :: IminS, ImaxS, JminS, JmaxS
142 integer,
intent(in) :: nrhs, nnew
146 real(r8),
intent(in) :: pmask(LBi:,LBj:)
148 real(r8),
intent(in) :: Hz(LBi:,LBj:,:)
149 real(r8),
intent(in) :: om_p(LBi:,LBj:)
150 real(r8),
intent(in) :: om_r(LBi:,LBj:)
151 real(r8),
intent(in) :: on_p(LBi:,LBj:)
152 real(r8),
intent(in) :: on_r(LBi:,LBj:)
153 real(r8),
intent(in) :: pm(LBi:,LBj:)
154 real(r8),
intent(in) :: pmon_p(LBi:,LBj:)
155 real(r8),
intent(in) :: pmon_r(LBi:,LBj:)
156 real(r8),
intent(in) :: pn(LBi:,LBj:)
157 real(r8),
intent(in) :: pnom_p(LBi:,LBj:)
158 real(r8),
intent(in) :: pnom_r(LBi:,LBj:)
159 real(r8),
intent(in) :: visc4_p(LBi:,LBj:)
160 real(r8),
intent(in) :: visc4_r(LBi:,LBj:)
162 real(r8),
intent(in) :: u(LBi:,LBj:,:,:)
163 real(r8),
intent(in) :: v(LBi:,LBj:,:,:)
172 real(r8),
intent(inout) :: ad_Hz(LBi:,LBj:,:)
173 real(r8),
intent(inout) :: ad_rufrc(LBi:,LBj:)
174 real(r8),
intent(inout) :: ad_rvfrc(LBi:,LBj:)
175 real(r8),
intent(inout) :: ad_u(LBi:,LBj:,:,:)
176 real(r8),
intent(inout) :: ad_v(LBi:,LBj:,:,:)
181 real(r8),
intent(in) :: pmask(LBi:UBi,LBj:UBj)
183 real(r8),
intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
184 real(r8),
intent(in) :: om_p(LBi:UBi,LBj:UBj)
185 real(r8),
intent(in) :: om_r(LBi:UBi,LBj:UBj)
186 real(r8),
intent(in) :: on_p(LBi:UBi,LBj:UBj)
187 real(r8),
intent(in) :: on_r(LBi:UBi,LBj:UBj)
188 real(r8),
intent(in) :: pm(LBi:UBi,LBj:UBj)
189 real(r8),
intent(in) :: pmon_p(LBi:UBi,LBj:UBj)
190 real(r8),
intent(in) :: pmon_r(LBi:UBi,LBj:UBj)
191 real(r8),
intent(in) :: pn(LBi:UBi,LBj:UBj)
192 real(r8),
intent(in) :: pnom_p(LBi:UBi,LBj:UBj)
193 real(r8),
intent(in) :: pnom_r(LBi:UBi,LBj:UBj)
194 real(r8),
intent(in) :: visc4_p(LBi:UBi,LBj:UBj)
195 real(r8),
intent(in) :: visc4_r(LBi:UBi,LBj:UBj)
197 real(r8),
intent(inout) :: u(LBi:UBi,LBj:UBj,N(ng),2)
198 real(r8),
intent(inout) :: v(LBi:UBi,LBj:UBj,N(ng),2)
207 real(r8),
intent(inout) :: ad_Hz(LBi:UBi,LBj:UBj,N(ng))
208 real(r8),
intent(inout) :: ad_rufrc(LBi:UBi,LBj:UBj)
209 real(r8),
intent(inout) :: ad_rvfrc(LBi:UBi,LBj:UBj)
210 real(r8),
intent(inout) :: ad_u(LBi:UBi,LBj:UBj,N(ng),2)
211 real(r8),
intent(inout) :: ad_v(LBi:UBi,LBj:UBj,N(ng),2)
216 integer :: IminU, IminV, ImaxU, ImaxV
217 integer :: JminU, JminV, JmaxU, JmaxV
220 real(r8) :: cff, cff1, cff2
221 real(r8) :: ad_cff, ad_cff1, ad_cff2
222 real(r8) :: adfac, adfac1, adfac2, adfac3, adfac4
224 real(r8),
dimension(IminS:ImaxS,JminS:JmaxS) :: LapU
225 real(r8),
dimension(IminS:ImaxS,JminS:JmaxS) :: LapV
226 real(r8),
dimension(IminS:ImaxS,JminS:JmaxS) :: UFe
227 real(r8),
dimension(IminS:ImaxS,JminS:JmaxS) :: VFe
228 real(r8),
dimension(IminS:ImaxS,JminS:JmaxS) :: UFx
229 real(r8),
dimension(IminS:ImaxS,JminS:JmaxS) :: VFx
231 real(r8),
dimension(IminS:ImaxS,JminS:JmaxS) :: ad_LapU
232 real(r8),
dimension(IminS:ImaxS,JminS:JmaxS) :: ad_LapV
233 real(r8),
dimension(IminS:ImaxS,JminS:JmaxS) :: ad_UFe
234 real(r8),
dimension(IminS:ImaxS,JminS:JmaxS) :: ad_VFe
235 real(r8),
dimension(IminS:ImaxS,JminS:JmaxS) :: ad_UFx
236 real(r8),
dimension(IminS:ImaxS,JminS:JmaxS) :: ad_VFx
238#include "set_bounds.h"
248 ad_lapu(imins:imaxs,jmins:jmaxs)=0.0_r8
249 ad_lapv(imins:imaxs,jmins:jmaxs)=0.0_r8
251 ad_ufe(imins:imaxs,jmins:jmaxs)=0.0_r8
252 ad_vfe(imins:imaxs,jmins:jmaxs)=0.0_r8
253 ad_ufx(imins:imaxs,jmins:jmaxs)=0.0_r8
254 ad_vfx(imins:imaxs,jmins:jmaxs)=0.0_r8
271 imaxu=min(iend+1,
lm(ng))
273 imaxv=min(iend+1,
lm(ng))
282 jmaxu=min(jend+1,
mm(ng))
284 jmaxv=min(jend+1,
mm(ng))
290 k_loop :
DO k=1,n(ng)
293 cff=visc4_r(i,j)*0.5_r8* &
295 & ((pn(i ,j)+pn(i+1,j))*u(i+1,j,k,nrhs)- &
296 & (pn(i-1,j)+pn(i ,j))*u(i ,j,k,nrhs))- &
298 & ((pm(i,j )+pm(i,j+1))*v(i,j+1,k,nrhs)- &
299 & (pm(i,j-1)+pm(i,j ))*v(i,j ,k,nrhs)))
300 ufx(i,j)=on_r(i,j)*on_r(i,j)*cff
301 vfe(i,j)=om_r(i,j)*om_r(i,j)*cff
306 cff=visc4_p(i,j)*0.5_r8* &
308 & ((pn(i ,j-1)+pn(i ,j))*v(i ,j,k,nrhs)- &
309 & (pn(i-1,j-1)+pn(i-1,j))*v(i-1,j,k,nrhs))+ &
311 & ((pm(i-1,j )+pm(i,j ))*u(i,j ,k,nrhs)- &
312 & (pm(i-1,j-1)+pm(i,j-1))*u(i,j-1,k,nrhs)))
316 ufe(i,j)=om_p(i,j)*om_p(i,j)*cff
317 vfx(i,j)=on_p(i,j)*on_p(i,j)*cff
325 lapu(i,j)=0.125_r8* &
326 & (pm(i-1,j)+pm(i,j))*(pn(i-1,j)+pn(i,j))* &
327 & ((pn(i-1,j)+pn(i,j))*(ufx(i,j )-ufx(i-1,j))+ &
328 & (pm(i-1,j)+pm(i,j))*(ufe(i,j+1)-ufe(i ,j)))
333 lapv(i,j)=0.125_r8* &
334 & (pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))* &
335 & ((pn(i,j-1)+pn(i,j))*(vfx(i+1,j)-vfx(i,j ))- &
336 & (pm(i,j-1)+pm(i,j))*(vfe(i ,j)-vfe(i,j-1)))
345 IF (
domain(ng)%Western_Edge(tile))
THEN
352 lapu(istr,j)=lapu(istr+1,j)
357 lapv(istr-1,j)=
gamma2(ng)*lapv(istr,j)
361 lapv(istr-1,j)=0.0_r8
368 IF (
domain(ng)%Eastern_Edge(tile))
THEN
371 lapu(iend+1,j)=0.0_r8
375 lapu(iend+1,j)=lapu(iend,j)
380 lapv(iend+1,j)=
gamma2(ng)*lapv(iend,j)
384 lapv(iend+1,j)=0.0_r8
391 IF (
domain(ng)%Southern_Edge(tile))
THEN
394 lapu(i,jstr-1)=
gamma2(ng)*lapu(i,jstr)
398 lapu(i,jstr-1)=0.0_r8
407 lapv(i,jstr)=lapv(i,jstr+1)
414 IF (
domain(ng)%Northern_Edge(tile))
THEN
417 lapu(i,jend+1)=
gamma2(ng)*lapu(i,jend)
421 lapu(i,jend+1)=0.0_r8
426 lapv(i,jend+1)=0.0_r8
430 lapv(i,jend+1)=lapv(i,jend)
438 IF (
domain(ng)%SouthWest_Corner(tile))
THEN
439 lapu(istr ,jstr-1)=0.5_r8* &
440 & (lapu(istr+1,jstr-1)+ &
442 lapv(istr-1,jstr )=0.5_r8* &
443 & (lapv(istr-1,jstr+1)+ &
450 IF (
domain(ng)%SouthEast_Corner(tile))
THEN
451 lapu(iend+1,jstr-1)=0.5_r8* &
452 & (lapu(iend ,jstr-1)+ &
453 & lapu(iend+1,jstr ))
454 lapv(iend+1,jstr )=0.5_r8* &
455 & (lapv(iend ,jstr )+ &
456 & lapv(iend+1,jstr+1))
462 IF (
domain(ng)%NorthWest_Corner(tile))
THEN
463 lapu(istr ,jend+1)=0.5_r8* &
464 & (lapu(istr+1,jend+1)+ &
466 lapv(istr-1,jend+1)=0.5_r8* &
467 & (lapv(istr ,jend+1)+ &
468 & lapv(istr-1,jend ))
474 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
475 lapu(iend+1,jend+1)=0.5_r8* &
476 & (lapu(iend ,jend+1)+ &
477 & lapu(iend+1,jend ))
478 lapv(iend+1,jend+1)=0.5_r8* &
479 & (lapv(iend ,jend+1)+ &
480 & lapv(iend+1,jend ))
489 cff=0.25_r8*(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
496 ad_cff2=ad_cff2-ad_v(i,j,k,nnew)
499 ad_cff1=ad_cff1-ad_rvfrc(i,j)
502 ad_cff1=ad_cff1+
dt(ng)*cff*ad_cff2
510 adfac1=adfac*(pn(i,j-1)+pn(i,j))
511 adfac2=adfac*(pm(i,j-1)+pm(i,j))
512 ad_vfx(i ,j )=ad_vfx(i ,j )-adfac1
513 ad_vfx(i+1,j )=ad_vfx(i+1,j )+adfac1
514 ad_vfe(i ,j-1)=ad_vfe(i ,j-1)+adfac2
515 ad_vfe(i ,j )=ad_vfe(i ,j )-adfac2
521 cff=0.25_r8*(pm(i-1,j)+pm(i,j))*(pn(i-1,j)+pn(i,j))
528 ad_cff2=ad_cff2-ad_u(i,j,k,nnew)
531 ad_cff1=ad_cff1-ad_rufrc(i,j)
534 ad_cff1=ad_cff1+
dt(ng)*cff*ad_cff2
542 adfac1=adfac*(pn(i-1,j)+pn(i,j))
543 adfac2=adfac*(pm(i-1,j)+pm(i,j))
544 ad_ufx(i-1,j )=ad_ufx(i-1,j )-adfac1
545 ad_ufx(i ,j )=ad_ufx(i ,j )+adfac1
546 ad_ufe(i ,j )=ad_ufe(i ,j )-adfac2
547 ad_ufe(i ,j+1)=ad_ufe(i ,j+1)+adfac2
559 ad_cff=ad_cff+on_p(i,j)*on_p(i,j)*ad_vfx(i,j)
563 ad_cff=ad_cff+om_p(i,j)*om_p(i,j)*ad_ufe(i,j)
568 ad_cff=ad_cff*pmask(i,j)
588 adfac=visc4_p(i,j)*0.125_r8*ad_cff
589 adfac1=adfac*(pmon_p(i,j)* &
590 & ((pn(i ,j-1)+pn(i ,j))*lapv(i ,j)- &
591 & (pn(i-1,j-1)+pn(i-1,j))*lapv(i-1,j))+ &
593 & ((pm(i-1,j )+pm(i,j ))*lapu(i,j )- &
594 & (pm(i-1,j-1)+pm(i,j-1))*lapu(i,j-1)))
595 adfac2=adfac*(hz(i-1,j ,k)+hz(i,j ,k)+ &
596 & hz(i-1,j-1,k)+hz(i,j-1,k))
597 adfac3=adfac2*pmon_p(i,j)
598 adfac4=adfac2*pnom_p(i,j)
599 ad_hz(i-1,j-1,k)=ad_hz(i-1,j-1,k)+adfac1
600 ad_hz(i-1,j ,k)=ad_hz(i-1,j ,k)+adfac1
601 ad_hz(i ,j-1,k)=ad_hz(i ,j-1,k)+adfac1
602 ad_hz(i ,j ,k)=ad_hz(i ,j ,k)+adfac1
603 ad_lapv(i-1,j)=ad_lapv(i-1,j)- &
604 & (pn(i-1,j-1)+pn(i-1,j))*adfac3
605 ad_lapv(i ,j)=ad_lapv(i ,j)+ &
606 & (pn(i ,j-1)+pn(i ,j))*adfac3
607 ad_lapu(i,j-1)=ad_lapu(i,j-1)- &
608 & (pm(i-1,j-1)+pm(i,j-1))*adfac4
609 ad_lapu(i,j )=ad_lapu(i,j )+ &
610 & (pm(i-1,j )+pm(i,j ))*adfac4
619 ad_cff=ad_cff+om_r(i,j)*om_r(i,j)*ad_vfe(i,j)
623 ad_cff=ad_cff+on_r(i,j)*on_r(i,j)*ad_ufx(i,j)
641 adfac=visc4_r(i,j)*0.5_r8*ad_cff
642 adfac1=adfac*hz(i,j,k)
643 adfac2=adfac1*pmon_r(i,j)
644 adfac3=adfac1*pnom_r(i,j)
645 ad_hz(i,j,k)=ad_hz(i,j,k)+ &
647 & ((pn(i ,j)+pn(i+1,j))*lapu(i+1,j)- &
648 & (pn(i-1,j)+pn(i ,j))*lapu(i ,j))- &
650 & ((pm(i,j )+pm(i,j+1))*lapv(i,j+1)- &
651 & (pm(i,j-1)+pm(i,j ))*lapv(i,j )))*adfac
652 ad_lapu(i ,j)=ad_lapu(i ,j)- &
653 & (pn(i-1,j)+pn(i ,j))*adfac2
654 ad_lapu(i+1,j)=ad_lapu(i+1,j)+ &
655 & (pn(i ,j)+pn(i+1,j))*adfac2
656 ad_lapv(i,j )=ad_lapv(i,j )+ &
657 & (pm(i,j-1)+pm(i,j ))*adfac3
658 ad_lapv(i,j+1)=ad_lapv(i,j+1)- &
659 & (pm(i,j )+pm(i,j+1))*adfac3
670 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
675 adfac=0.5_r8*ad_lapv(iend+1,jend+1)
676 ad_lapv(iend ,jend+1)=ad_lapv(iend ,jend+1)+adfac
677 ad_lapv(iend+1,jend )=ad_lapv(iend+1,jend )+adfac
678 ad_lapv(iend+1,jend+1)=0.0_r8
683 adfac=0.5_r8*ad_lapu(iend+1,jend+1)
684 ad_lapu(iend ,jend+1)=ad_lapu(iend ,jend+1)+adfac
685 ad_lapu(iend+1,jend )=ad_lapu(iend+1,jend )+adfac
686 ad_lapu(iend+1,jend+1)=0.0_r8
692 IF (
domain(ng)%NorthWest_Corner(tile))
THEN
697 adfac=0.5_r8*ad_lapv(istr-1,jend+1)
698 ad_lapv(istr ,jend+1)=ad_lapv(istr ,jend+1)+adfac
699 ad_lapv(istr-1,jend )=ad_lapv(istr-1,jend )+adfac
700 ad_lapv(istr-1,jend+1)=0.0_r8
705 adfac=0.5_r8*ad_lapu(istr,jend+1)
706 ad_lapu(istr+1,jend+1)=ad_lapu(istr+1,jend+1)+adfac
707 ad_lapu(istr ,jend )=ad_lapu(istr ,jend )+adfac
708 ad_lapu(istr ,jend+1)=0.0_r8
714 IF (
domain(ng)%SouthEast_Corner(tile))
THEN
719 adfac=0.5_r8*ad_lapv(iend+1,jstr)
720 ad_lapv(iend ,jstr )=ad_lapv(iend ,jstr )+adfac
721 ad_lapv(iend+1,jstr+1)=ad_lapv(iend+1,jstr+1)+adfac
722 ad_lapv(iend+1,jstr )=0.0_r8
727 adfac=0.5_r8*ad_lapu(iend+1,jstr-1)
728 ad_lapu(iend ,jstr-1)=ad_lapu(iend ,jstr-1)+adfac
729 ad_lapu(iend+1,jstr )=ad_lapu(iend+1,jstr )+adfac
730 ad_lapu(iend+1,jstr-1)=0.0_r8
736 IF (
domain(ng)%SouthWest_Corner(tile))
THEN
741 adfac=0.5_r8*ad_lapv(istr-1,jstr)
742 ad_lapv(istr-1,jstr+1)=ad_lapv(istr-1,jstr+1)+adfac
743 ad_lapv(istr ,jstr )=ad_lapv(istr ,jstr )+adfac
744 ad_lapv(istr-1,jstr )=0.0_r8
749 adfac=0.5_r8*ad_lapu(istr,jstr-1)
750 ad_lapu(istr+1,jstr-1)=ad_lapu(istr+1,jstr-1)+adfac
751 ad_lapu(istr ,jstr )=ad_lapu(istr ,jstr )+adfac
752 ad_lapu(istr ,jstr-1)=0.0_r8
757 IF (
domain(ng)%Northern_Edge(tile))
THEN
762 ad_lapv(i,jend+1)=0.0_r8
768 ad_lapv(i,jend)=ad_lapv(i,jend)+ad_lapv(i,jend+1)
769 ad_lapv(i,jend+1)=0.0_r8
776 ad_lapu(i,jend)=ad_lapu(i,jend)+ &
777 &
gamma2(ng)*ad_lapu(i,jend+1)
778 ad_lapu(i,jend+1)=0.0_r8
784 ad_lapu(i,jend+1)=0.0_r8
791 IF (
domain(ng)%Southern_Edge(tile))
THEN
796 ad_lapv(i,jstrv-1)=0.0_r8
802 ad_lapv(i,jstrv)=ad_lapv(i,jstrv)+ad_lapv(i,jstrv-1)
803 ad_lapv(i,jstrv-1)=0.0_r8
810 ad_lapu(i,jstr)=ad_lapu(i,jstr)+ &
811 &
gamma2(ng)*ad_lapu(i,jstr-1)
812 ad_lapu(i,jstr-1)=0.0_r8
818 ad_lapu(i,jstr-1)=0.0_r8
825 IF (
domain(ng)%Eastern_Edge(tile))
THEN
830 ad_lapv(iend,j)=ad_lapv(iend,j)+ &
831 &
gamma2(ng)*ad_lapv(iend+1,j)
832 ad_lapv(iend+1,j)=0.0_r8
838 ad_lapv(iend+1,j)=0.0_r8
845 ad_lapu(iend+1,j)=0.0_r8
851 ad_lapu(iend,j)=ad_lapu(iend,j)+ad_lapu(iend+1,j)
852 ad_lapu(iend+1,j)=0.0_r8
859 IF (
domain(ng)%Western_Edge(tile))
THEN
864 ad_lapv(istr,j)=ad_lapv(istr,j)+ &
865 &
gamma2(ng)*ad_lapv(istr-1,j)
866 ad_lapv(istr-1,j)=0.0_r8
872 ad_lapv(istr-1,j)=0.0_r8
879 ad_lapu(istru-1,j)=0.0_r8
885 ad_lapu(istru,j)=ad_lapu(istru,j)+ad_lapu(istru-1,j)
886 ad_lapu(istru-1,j)=0.0_r8
896 cff=0.125_r8*(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
903 adfac=cff*ad_lapv(i,j)
904 adfac1=adfac*(pn(i,j-1)+pn(i,j))
905 adfac2=adfac*(pm(i,j-1)+pm(i,j))
906 ad_vfx(i ,j)=ad_vfx(i ,j)-adfac1
907 ad_vfx(i+1,j)=ad_vfx(i+1,j)+adfac1
908 ad_vfe(i,j-1)=ad_vfe(i,j-1)+adfac2
909 ad_vfe(i, j)=ad_vfe(i, j)-adfac2
915 cff=0.125_r8*(pm(i-1,j)+pm(i,j))*(pn(i-1,j)+pn(i,j))
922 adfac=cff*ad_lapu(i,j)
923 adfac1=adfac*(pn(i-1,j)+pn(i,j))
924 adfac2=adfac*(pm(i-1,j)+pm(i,j))
925 ad_ufx(i-1,j)=ad_ufx(i-1,j)-adfac1
926 ad_ufx(i ,j )=ad_ufx(i ,j)+adfac1
927 ad_ufe(i,j )=ad_ufe(i,j )-adfac2
928 ad_ufe(i,j+1)=ad_ufe(i,j+1)+adfac2
944 ad_cff=ad_cff+on_p(i,j)*on_p(i,j)*ad_vfx(i,j)
948 ad_cff=ad_cff+om_p(i,j)*om_p(i,j)*ad_ufe(i,j)
953 ad_cff=ad_cff*pmask(i,j)
963 adfac=visc4_p(i,j)*0.5_r8*ad_cff
964 adfac1=adfac*pmon_p(i,j)
965 adfac2=adfac*pnom_p(i,j)
966 ad_v(i-1,j,k,nrhs)=ad_v(i-1,j,k,nrhs)- &
967 & (pn(i-1,j-1)+pn(i-1,j))*adfac1
968 ad_v(i ,j,k,nrhs)=ad_v(i ,j,k,nrhs)+ &
969 & (pn(i ,j-1)+pn(i ,j))*adfac1
970 ad_u(i,j ,k,nrhs)=ad_u(i,j ,k,nrhs)+ &
971 & (pm(i-1,j )+pm(i,j ))*adfac2
972 ad_u(i,j-1,k,nrhs)=ad_u(i,j-1,k,nrhs)- &
973 & (pm(i-1,j-1)+pm(i,j-1))*adfac2
981 ad_cff=ad_cff+om_r(i,j)*om_r(i,j)*ad_vfe(i,j)
985 ad_cff=ad_cff+on_r(i,j)*on_r(i,j)*ad_ufx(i,j)
995 adfac=visc4_r(i,j)*0.5_r8*ad_cff
996 adfac1=adfac*pmon_r(i,j)
997 adfac2=adfac*pnom_r(i,j)
998 ad_u(i ,j,k,nrhs)=ad_u(i ,j,k,nrhs)- &
999 & (pn(i-1,j)+pn(i ,j))*adfac1
1000 ad_u(i+1,j,k,nrhs)=ad_u(i+1,j,k,nrhs)+ &
1001 & (pn(i ,j)+pn(i+1,j))*adfac1
1002 ad_v(i,j ,k,nrhs)=ad_v(i,j ,k,nrhs)+ &
1003 & (pm(i,j-1)+pm(i,j ))*adfac2
1004 ad_v(i,j+1,k,nrhs)=ad_v(i,j+1,k,nrhs)- &
1005 & (pm(i,j )+pm(i,j+1))*adfac2