SUBROUTINE PGCOR C ****************************************************************** C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . C SUBPROGRAM: PGCOR PRESSURE GRADIENT/CORIOLIS CALC C PRGRMMR: JANJIC ORG: W/NP22 DATE: 93-10-28 C C ABSTRACT: C PGCOR CALCULATES THE PRESSURE GRADIENT FORCE, UPDATES THE C VELOCITY COMPONENTS DUE TO THE EFFECT OF THE PRESSURE GRADIENT C AND CORIOLIS FORCES. C C PROGRAM HISTORY LOG: C 87-06-?? JANJIC - ORIGINATOR C 95-03-25 BLACK - CONVERSION FROM 1-D TO 2-D IN HORIZONTAL C 96-03-29 BLACK - ADDED EXTERNAL EDGE C 97-03-17 MESINGER - SPLIT FROM PFDHT C 98-10-28 BLACK - MODIFIED FOR DISTRIBUTED MEMORY C C USAGE: CALL PFDHT FROM MAIN PROGRAM EBU C INPUT ARGUMENT LIST: C NONE C C OUTPUT ARGUMENT LIST: C NONE C C OUTPUT FILES: C NONE C C SUBPROGRAMS CALLED: C C UNIQUE: NONE C C LIBRARY: NONE C C COMMON BLOCKS: CTLBLK C MASKS C DYNAM C VRBLS C CONTIN Cligia NHYDRO C INDX C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE : IBM SP C$$$ C*********************************************************************** INCLUDE "parmeta" INCLUDE "mpp.h" #include "sp.h" C----------------------------------------------------------------------- P A R A M E T E R & (IMJM=IM*JM-JM/2,JAM=6+2*(JM-10),LP1=LM+1) C----------------------------------------------------------------------- Cligia nhydro L O G I C A L & RUN,FIRST,RESTRT,SIGMA C---------------------------------------------------------------------- INCLUDE "CTLBLK.comm" C----------------------------------------------------------------------- INCLUDE "MASKS.comm" C----------------------------------------------------------------------- INCLUDE "DYNAM.comm" C----------------------------------------------------------------------- INCLUDE "VRBLS.comm" C----------------------------------------------------------------------- INCLUDE "PVRBLS.comm" C----------------------------------------------------------------------- INCLUDE "CLDWTR.comm" C----------------------------------------------------------------------- INCLUDE "CONTIN.comm" c----------------------------------------------------------------------- INCLUDE "NHYDRO.comm" C----------------------------------------------------------------------- INCLUDE "INDX.comm" C----------------------------------------------------------------------- Cligia R E A L & PINTLG(IDIM1:IDIM2,JDIM1:JDIM2,LM+1) Cligia D I M E N S I O N & FILO (IDIM1:IDIM2,JDIM1:JDIM2),RDPD (IDIM1:IDIM2,JDIM1:JDIM2) &,ADPDX (IDIM1:IDIM2,JDIM1:JDIM2),RDPDX (IDIM1:IDIM2,JDIM1:JDIM2) &,ADPDY (IDIM1:IDIM2,JDIM1:JDIM2),RDPDY (IDIM1:IDIM2,JDIM1:JDIM2) &,FIUP (IDIM1:IDIM2,JDIM1:JDIM2),ADPDNE(IDIM1:IDIM2,JDIM1:JDIM2) &,ADPDSE(IDIM1:IDIM2,JDIM1:JDIM2) &,PEW (IDIM1:IDIM2,JDIM1:JDIM2),PNS (IDIM1:IDIM2,JDIM1:JDIM2) &,PCEW (IDIM1:IDIM2,JDIM1:JDIM2),PCNS (IDIM1:IDIM2,JDIM1:JDIM2) &,DPFEW (IDIM1:IDIM2,JDIM1:JDIM2),DPFNS (IDIM1:IDIM2,JDIM1:JDIM2) &,FNS (IDIM1:IDIM2,JDIM1:JDIM2),TNS (IDIM1:IDIM2,JDIM1:JDIM2) &,DPNE (IDIM1:IDIM2,JDIM1:JDIM2),DPSE (IDIM1:IDIM2,JDIM1:JDIM2) &,DCNE (IDIM1:IDIM2,JDIM1:JDIM2),DCSE (IDIM1:IDIM2,JDIM1:JDIM2) &,DPFNE (IDIM1:IDIM2,JDIM1:JDIM2),DPFSE (IDIM1:IDIM2,JDIM1:JDIM2) &,VM (IDIM1:IDIM2,JDIM1:JDIM2),F0 (IDIM1:IDIM2,JDIM1:JDIM2) &,UP (IDIM1:IDIM2,JDIM1:JDIM2),VP (IDIM1:IDIM2,JDIM1:JDIM2) &,PVNE (IDIM1:IDIM2,JDIM1:JDIM2),PVSE (IDIM1:IDIM2,JDIM1:JDIM2) &,HM (IDIM1:IDIM2,JDIM1:JDIM2) C D I M E N S I O N & DPDE(IDIM1:IDIM2,JDIM1:JDIM2),FIM (IDIM1:IDIM2,JDIM1:JDIM2) &,APEL(IDIM1:IDIM2,JDIM1:JDIM2),PCXC(IDIM1:IDIM2,JDIM1:JDIM2) &,UDY (IDIM1:IDIM2,JDIM1:JDIM2),VDX (IDIM1:IDIM2,JDIM1:JDIM2) &,TEW (IDIM1:IDIM2,JDIM1:JDIM2),FEW (IDIM1:IDIM2,JDIM1:JDIM2) &,TNE (IDIM1:IDIM2,JDIM1:JDIM2),TSE (IDIM1:IDIM2,JDIM1:JDIM2) &,FNE (IDIM1:IDIM2,JDIM1:JDIM2),FSE (IDIM1:IDIM2,JDIM1:JDIM2) &,PNE (IDIM1:IDIM2,JDIM1:JDIM2),PSE (IDIM1:IDIM2,JDIM1:JDIM2) &,CNE (IDIM1:IDIM2,JDIM1:JDIM2),CSE (IDIM1:IDIM2,JDIM1:JDIM2) &,PPNE(IDIM1:IDIM2,JDIM1:JDIM2),PPSE(IDIM1:IDIM2,JDIM1:JDIM2) &,PCNE(IDIM1:IDIM2,JDIM1:JDIM2),PCSE(IDIM1:IDIM2,JDIM1:JDIM2) Cligia &,ALP1 (IDIM1:IDIM2,JDIM1:JDIM2) &,DFDZ (IDIM1:IDIM2,JDIM1:JDIM2) Cligia C----------------------------------------------------------------------- C----------------------------------------------------------------------- Cligia ALP1, ADPDX, ADPDY, DFDZ CALL ZERO2(ALP1) CALL ZERO2(DPDE) CALL ZERO2(FIM) CALL ZERO2(APEL) CALL ZERO2(ADPDX) CALL ZERO2(ADPDY) CALL ZERO2(DFDZ) CALL ZERO2(PCXC) CALL ZERO2(UDY) CALL ZERO2(VDX) CALL ZERO2(TEW) CALL ZERO2(FEW) CALL ZERO2(TNE) CALL ZERO2(TSE) CALL ZERO2(FNE) CALL ZERO2(FSE) CALL ZERO2(PNE) CALL ZERO2(CNE) CALL ZERO2(PSE) CALL ZERO2(CSE) CALL ZERO2(PPNE) CALL ZERO2(PCNE) CALL ZERO2(PPSE) CALL ZERO2(PCSE) call zero2(pdsl) call zero2(adpdne) C----------------------------------------------------------------------- C--------------PREPARATORY CALCULATIONS--------------------------------- C----------------------------------------------------------------------- Cligia IF(SIGMA)THEN !$omp parallel do DO 50 J=MYJS_P5,MYJE_P5 DO 50 I=MYIS_P5,MYIE_P5 FILO(I,J)=FIS(I,J) PDSL(I,J)=PD(I,J) 50 CONTINUE ELSE !$omp parallel do DO 100 J=MYJS_P5,MYJE_P5 DO 100 I=MYIS_P5,MYIE_P5 FILO(I,J)=0.0 PDSL(I,J)=RES(I,J)*PD(I,J) 100 CONTINUE ENDIF C IF(HYDRO)THEN !$omp parallel do DO L=1,LM+1 DO J=MYJS_P5,MYJE_P5 DO I=MYIS_P5,MYIE_P5 PINTLG(I,J,L)=ALOG(ETA(L)*PDSL(I,J)+PT) ENDDO ENDDO ENDDO ELSE !$omp parallel do DO L=1,LM+1 DO J=MYJS_P5,MYJE_P5 DO I=MYIS_P5,MYIE_P5 PINTLG(I,J,L)=ALOG(PINT(I,J,L)) if (PINTLG(I,J,L) .eq. PINTLG(I,J,L-1)) then write(6,*) 'same pintlg at different levels: ', MYPE,I,J,L STOP 690 write(6,*) 'PINT vals: ', PINT(I,J,L),PINT(I,J,L-1), + PINT(I,J,L-2),PINT(I,J,L-3) write(6,*) 'PINT vals at I-1: ', PINT(I-1,J,L),PINT(I-1,J,L-1), + PINT(I-1,J,L-2),PINT(I-1,J,L-3) endif ENDDO ENDDO ENDDO ENDIF Cligia C----------------------------------------------------------------------- Cligia !$omp parallel do private (alp1p) DO J=MYJS_P5,MYJE_P5 DO I=MYIS_P5,MYIE_P5 ALP1P=PINTLG(I,J,LM+1) ALP1(I,J)=ALP1P ENDDO ENDDO Cligia C !$omp parallel do DO J=MYJS_P4,MYJE_P4 DO I=MYIS_P4,MYIE_P4 ADPDX(I,J)=0. ADPDY(I,J)=0. ENDDO ENDDO C--------------MAIN VERTICAL INTEGRATION LOOP--------------------------- cligia FIM=0 cligia DO 400 L=LM,1,-1 C----------------------------------------------------------------------- Cligia C*** C*** INTEGRATE THE GEOPOTENTIAL C*** !$omp parallel do private (alp1p,dfi,fiupk,rdpds) DO 125 J=MYJS_P5,MYJE_P5 DO 125 I=MYIS_P5,MYIE_P5 C ALP1P=PINTLG(I,J,L) C Cdule DFI=(Q(I,J,L)*0.608+1.)*T(I,J,L)*R*(ALP1(I,J)-ALP1P)/DWDT(I,J,L) DFI=(Q(I,J,L)*0.608+1.)*T(I,J,L)*R*(ALP1(I,J)-ALP1P) & /(1+CWM(I,J,L))/DWDT(I,J,L) if (abs(DFI) .lt. 2.e13) then else write(6,*) 'BAD DFI: ', DFI write(6,*) 'Q,T: ', Q(I,J,L),T(I,J,L) write(6,*) 'ALP vals: ', ALP1(I,J),ALP1P write(6,*) 'DWDT= ', DWDT(I,J,L) endif C RDPDS=1./(DETA(L)*PDSL(I,J)) FIUPK=FILO(I,J)+DFI if (abs(FIUPK) .lt. 2.e13) then else write(6,*) 'bad FIUPK. FILO, DFI ', FILO(I,J), DFI endif FIM(I,J)=FILO(I,J)+FIUPK C FILO(I,J)=(FIUPK-DFL(L))*HTM(I,J,L)+DFL(L) C if (MYPE.eq.4.and.I.eq.1.and.J.eq.2.and.L.eq.25) then C write(6,*) 'L,FILO, DFI-->: ', L,FILO(I,J),DFI C write(6,*) 'T,Q,ALP1,ALP1P,DWDT: ', T(I,J,L),Q(I,J,L),ALP1(I,J), C + ALP1P,DWDT(I,J,L) C endif if (abs(FILO(I,J)) .lt. 20000000.) then else write(6,*) 'bad FILO value ', FILO(I,J),' on PE: ' , + MYPE, 'at ', I,J write(6,*) 'FIUPK,DFL: ', FIUPK, DFL(L) STOP 999 endif ALP1(I,J)=ALP1P 125 CONTINUE Cligia C----------------------------------------------------------------------- Cligia !$omp parallel do private (alp1p,alp1pl,alp2p,alp2pl,dfi) DO 205 J=MYJS_P5,MYJE_P5 DO 205 I=MYIS_P5,MYIE_P5 HM(I,J)=HTM(I,J,L)*HBM2(I,J) VM(I,J)=VTM(I,J,L)*VBM2(I,J) C ALP1P =PINTLG(I,J,L) ALP1PL=PINTLG(I,J,L+1) ALP2P =ALP1P*ALP1P ALP2PL=ALP1PL*ALP1PL C Cdule DFI=(Q(I,J,L)*0.608+1.)*T(I,J,L)*R*(ALP1PL-ALP1P)/DWDT(I,J,L) DFI=(Q(I,J,L)*0.608+1.)*T(I,J,L)*R*(ALP1PL-ALP1P) & /(1+CWM(I,J,L))/DWDT(I,J,L) if (abs(DFI) .le. 2.e13) then else write(6,*) 'BAD DFI:::: ', DFI endif if (abs(DWDT(I,J,L)) .le. 2.e13) then else write(6,*) 'BAD DWDTI:::: ', DWDT(I,J,L) endif if (abs(ALP2PL) .le. 2.e13) then else write(6,*) 'BAD ALP2PL:::: ', ALP2PL endif if (abs(ALP2P) .le. 2.e13) then else write(6,*) 'BAD ALP2P:::: ', ALP2P endif DFDZ(I,J)=DFI*DWDT(I,J,L)/(ALP2PL-ALP2P) if (abs(DFDZ(I,J)) .le. 2.e13) then else write(6,*) 'on PE: ', MYPE write(6,*) 'at = ', I,J,L write(6,*) 'DFDZ= ', DFDZ(I,J) write(6,*) 'DFI= ', DFI write(6,*) 'DWDT= ', DWDT(I,J,L) write(6,*) 'denom= ', ALP2PL-ALP2P write(6,*) 'PINTLG(L) ', PINTLG(I,J,L) write(6,*) 'PINTLG(L+1) ', PINTLG(I,J,L+1) endif APEL(I,J)=(ALP2PL+ALP2P)*0.5 DPDE(I,J)=DETA(L)*PDSL(I,J) 205 CONTINUE Cligia C !$omp parallel do cligia DO 215 J=MYJS,MYJE cligia DO 215 I=MYIS,MYIE DO 215 J=MYJS_P1,MYJE_P1 DO 215 I=MYIS_P1,MYIE_P1 RDPD(I,J)=1./DPDE(I,J) 215 CONTINUE C !$omp parallel do cligia DO 220 J=MYJS1_P2,MYJE1_P2 cligia DO 220 I=MYIS_P2,MYIE1_P2 DO 220 J=MYJS1_P3,MYJE1_P3 DO 220 I=MYIS_P3,MYIE_P3 ADPDX(I,J)=DPDE(I+IVW(J),J)+DPDE(I+IVE(J),J) ADPDY(I,J)=DPDE(I,J-1)+DPDE(I,J+1) RDPDX(I,J)=1./ADPDX(I,J) RDPDY(I,J)=1./ADPDY(I,J) 220 CONTINUE C----------------------------------------------------------------------- Cligia !$omp parallel do c DO 230 J=MYJS_P4,MYJE_P4 c DO 230 I=MYIS_P4,MYIE_P4 c APEL(I,J)=PT+AETA(L)*PDSL(I,J) c RTOP(I,J,L)=R*T(I,J,L)*(1.+0.608*Q(I,J,L))/APEL(I,J) c FIUPK=FILO(I,J)+RTOP(I,J,L)*DPDE(I,J) c FIM(I,J)=FILO(I,J)+FIUPK c FILO(I,J)=DFL(L)+HTM(I,J,L)*(FIUPK-DFL(L)) c 230 CONTINUE Cligia C--------------DIAGONAL CONTRIBUTIONS TO PRESSURE GRADIENT FORCE-------- C Cligia CC Having problem with global boundary here CC !$omp parallel do if (MYPE .eq. 5 .and. L .eq. LM ) then I=MYIE_P4 J=MYJE_P4-2 C write(6,*) 'bad FIM point: ',I+IHE(J),J+1,FIM (I+IHE(J),J+1) endif DO 240 J=MYJS_P4,MYJE_P4 DO 240 I=MYIS_P4,MYIE_P4 ADPDNE(I,J)=DPDE(I+IHE(J),J+1)+DPDE(I,J) if ( abs(FIM (I+IHE(J),J+1)) .lt. 2000000.) then else write(6,*) 'using FIM val: ', FIM (I+IHE(J),J+1), 'at point', + I+IHE(J),J+1, 'on PE: ', MYPE endif PNE(I,J)=(FIM (I+IHE(J),J+1)-FIM (I,J)) 1 *(DWDT(I+IHE(J),J+1,L)+DWDT(I,J,L)) PPNE(I,J)=PNE(I,J)*ADPDNE(I,J) CNE(I,J)=(DFDZ(I+IHE(J),J+1)+DFDZ(I,J))*2. 1 *(APEL(I+IHE(J),J+1)-APEL(I,J)) PCNE(I,J)=CNE(I,J)*ADPDNE(I,J) 240 CONTINUE Cligia !$omp parallel do Cligia DO 250 J=MYJS1_P3,MYJE_P3 Cligia DO 250 I=MYIS_P3,MYIE1_P3 DO 250 J=MYJS1_P4,MYJE_P4 DO 250 I=MYIS_P4,MYIE1_P4 ADPDSE(I,J)=DPDE(I+IHE(J),J-1)+DPDE(I,J) cligia PSE(I,J)=2.*(FIM(I+IHE(J),J-1)-FIM(I,J)) cligia PSE(I,J)=(FIM (I+IHE(J),J-1)-FIM (I,J)) 1 *(DWDT(I+IHE(J),J-1,L)+DWDT(I,J,L)) cligia PPSE(I,J)=PSE(I,J)*ADPDSE(I,J) Cligia CSE(I,J)=2.*(RTOP(I+IHE(J),J-1,L)+RTOP(I,J,L)) Cligia 1 *(APEL(I+IHE(J),J-1)-APEL(I,J)) CSE(I,J)=(DFDZ(I+IHE(J),J-1)+DFDZ(I,J))*2. 1 *(APEL(I+IHE(J),J-1)-APEL(I,J)) PCSE(I,J)=CSE(I,J)*ADPDSE(I,J) 250 CONTINUE C--------------LAT & LONG PRESSURE FORCE COMPONENTS--------------------- !$omp parallel do private(dcnek,dcsek,dpnek,dpsek) Cligia DO 280 J=MYJS1_P2,MYJE1_P2 Cligia DO 280 I=MYIS_P2,MYIE_P2 DO 280 J=MYJS1_P3,MYJE1_P3 DO 280 I=MYIS_P3,MYIE_P3 DPNEK=PNE(I+IVW(J),J)+PNE(I,J-1) DPSEK=PSE(I+IVW(J),J)+PSE(I,J+1) PEW(I,J)=DPNEK+DPSEK PNS(I,J)=DPNEK-DPSEK DCNEK=CNE(I+IVW(J),J)+CNE(I,J-1) DCSEK=CSE(I+IVW(J),J)+CSE(I,J+1) PCEW(I,J)=(DCNEK+DCSEK)*ADPDX(I,J) PCNS(I,J)=(DCNEK-DCSEK)*ADPDY(I,J) 280 CONTINUE C--------------UPDATE U AND V (CORIOLIS & PGF)-------------------------- !$omp parallel do private(dpfnek,dpfsek) Cligia DO 290 J=MYJS2_P2,MYJE2_P2 Cligia DO 290 I=MYIS_P2,MYIE1_P2 DO 290 J=MYJS2_P3,MYJE2_P3 DO 290 I=MYIS_P3,MYIE1_P3 DPFNEK=((PPNE(I+IVW(J),J)+PPNE(I,J-1)) 1 +(PCNE(I+IVW(J),J)+PCNE(I,J-1)))*2. DPFSEK=((PPSE(I+IVW(J),J)+PPSE(I,J+1)) 1 +(PCSE(I+IVW(J),J)+PCSE(I,J+1)))*2. DPFEW(I,J)=DPFNEK+DPFSEK DPFNS(I,J)=DPFNEK-DPFSEK 290 CONTINUE C !$omp parallel do private(f0k,upk,utk,vpk,vtk) DO 300 J=MYJS2_P2,MYJE2_P2 DO 300 I=MYIS_P2,MYIE1_P2 F0K=U(I,J,L)*CURV(I,J)+F(I,J) VM(I,J)=VTM(I,J,L)*VBM2(I,J) UPK=((DPFEW(I,J)+PCEW(I,J))*RDPDX(I,J) 1 +PEW(I,J))*CPGFU(I,J)+F0K*V(I,J,L)+U(I,J,L) VPK=((DPFNS(I,J)+PCNS(I,J))*RDPDY(I,J) 1 +PNS(I,J))*CPGFV-F0K*U(I,J,L)+V(I,J,L) UTK=U(I,J,L) VTK=V(I,J,L) U(I,J,L)=((F0K*VPK+UPK)/(F0K*F0K+1.)-U(I,J,L)) 1 *VM(I,J)+U(I,J,L) V(I,J,L)=(VPK-F0K*U(I,J,L)-V(I,J,L)) 1 *VM(I,J)+V(I,J,L) 300 CONTINUE C----------------------------------------------------------------------- 400 CONTINUE C----------------------------------------------------------------------- RETURN END