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 00-10-20 BLACK - INCORPORATED PRESSURE GRADIENT METHOD C FROM MESO MODEL C C USAGE: CALL PGCOR 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 LOOPS C DYNAM C VRBLS C CONTIN C NHYDRO C INDX C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE : IBM SP C$$$ C*********************************************************************** C----------------------------------------------------------------------- INCLUDE "parmeta" INCLUDE "mpp.h" C----------------------------------------------------------------------- P A R A M E T E R & (LP1=LM+1,JAM=6+2*(JM-10)) C----------------------------------------------------------------------- L O G I C A L & RUN,FIRST,RESTRT,SIGMA C---------------------------------------------------------------------- INCLUDE "CTLBLK.comm" c----------------------------------------------------------------------- include "LOOPS.comm" C----------------------------------------------------------------------- INCLUDE "MASKS.comm" C----------------------------------------------------------------------- INCLUDE "INDX.comm" c----------------------------------------------------------------------- INCLUDE "DYNAM.comm" C----------------------------------------------------------------------- INCLUDE "VRBLS.comm" c----------------------------------------------------------------------- INCLUDE "CONTIN.comm" c----------------------------------------------------------------------- INCLUDE "NHYDRO.comm" c----------------------------------------------------------------------- INCLUDE "CLDWTR.comm" c----------------------------------------------------------------------- REAL & PINTLG(IDIM1:IDIM2,JDIM1:JDIM2,LM+1) C C REAL & FIM (IDIM1:IDIM2,JDIM1:JDIM2) &,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) &,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) &,HM (IDIM1:IDIM2,JDIM1:JDIM2),VM (IDIM1:IDIM2,JDIM1:JDIM2) C REAL & DPDE (IDIM1:IDIM2,JDIM1:JDIM2) &,APEL (IDIM1:IDIM2,JDIM1:JDIM2) &,ALP1 (IDIM1:IDIM2,JDIM1:JDIM2) &,DFDZ (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) C----------------------------------------------------------------------- C----------------------------------------------------------------------- CHOU REAL SPINT, RIS CALL ZERO2(ALP1) CALL ZERO2(DPDE) CALL ZERO2(APEL) CALL ZERO2(ADPDX) CALL ZERO2(ADPDY) CALL ZERO2(DFDZ) CALL ZERO2(PNE) CALL ZERO2(PSE) CALL ZERO2(CNE) CALL ZERO2(CSE) CALL ZERO2(PPNE) CALL ZERO2(PPSE) CALL ZERO2(PCNE) CALL ZERO2(PCSE) C----------------------------------------------------------------------- C--------------PREPARATORY CALCULATIONS--------------------------------- C----------------------------------------------------------------------- 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,*) 'same1 pintlg at different levels: ', MYPE,I,J,L C write(6,*) 'PINT1 vals: ', PINT(I,J,L),PINT(I,J,L-1), + PINT(I,J,L-2),PINT(I,J,L-3) C write(6,*) 'PINT1 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) CHOU if equal PINTs replace by 4 point average write(6,*) '>>PINT1 vals at I-1, I+1, J-1, J+1 (L-1):', + PINT(I-1,J,L-1),PINT(I+1,J,L-1),PINT(I,J-1,L-1),PINT(I,J+1,L-1) SPINT=(PINT(I-1,J,L-1)+PINT(I+1,J,L-1)+PINT(I,J-1,L-1)+ + PINT(I,J+1,L-1)) is=0 if (PINT(I-1,J,L-1).ge.2000.0) is=is+1 if (PINT(I+1,J,L-1).ge.2000.0) is=is+1 if (PINT(I,J-1,L-1).ge.2000.0) is=is+1 if (PINT(I,J+1,L-1).ge.2000.0) is=is+1 ris=1./is PINT(I,J,L-1)= SPINT*ris PINTLG(I,J,L-1)=ALOG(PINT(I,J,L-1)) write(6,*) '>>mype1,i,j,l,PINT(i,j,l):',mype,i,j,l, + PINT(i,j,l),PINT(i,j,l-1),PINTLG(i,j,l),PINTLG(i,j,l-1) CGUSTAVO if (PINTLG(I,J,L) .eq. PINTLG(I,J,L-1)) then CHOU if equal PINTs replace by 4 point average write(6,*) '>>PINT1 vals1 at I-1, I+1, J-1, J+1 (L-1):', + PINT(I-1,J,L-1),PINT(I+1,J,L-1),PINT(I,J-1,L-1),PINT(I,J+1,L-1), + PINT(I-1,J-1,L-1),PINT(I-1,J+1,L-1),PINT(I+1,J-1,L-1), + PINT(I+1,J+1,L-1) SPINT=(PINT(I-1,J,L-1)+PINT(I+1,J,L-1)+PINT(I,J-1,L-1)+ + PINT(I,J+1,L-1)+PINT(I-1,J-1,L-1)+PINT(I-1,J+1,L-1)+ + PINT(I+1,J-1,L-1)+PINT(I+1,J+1,L-1)) is=0 if (PINT(I-1,J,L-1).ge.2000.0) is=is+1 if (PINT(I+1,J,L-1).ge.2000.0) is=is+1 if (PINT(I,J-1,L-1).ge.2000.0) is=is+1 if (PINT(I,J+1,L-1).ge.2000.0) is=is+1 if (PINT(I-1,J-1,L-1).ge.2000.0) is=is+1 if (PINT(I-1,J+1,L-1).ge.2000.0) is=is+1 if (PINT(I+1,J-1,L-1).ge.2000.0) is=is+1 if (PINT(I+1,J+1,L-1).ge.2000.0) is=is+1 ris=1./is PINT(I,J,L-1)= SPINT*ris PINTLG(I,J,L-1)=ALOG(PINT(I,J,L-1)) write(6,*) '>>mype1.1,i,j,l,PINT(i,j,l):',mype,i,j,l, + PINT(i,j,l),PINT(i,j,l-1),PINTLG(i,j,l),PINTLG(i,j,l-1) endif CHOU STOP 690 endif ENDDO ENDDO ENDDO ENDIF C----------------------------------------------------------------------- !$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 C----------------------------------------------------------------------- C-------------------- MAIN VERTICAL INTEGRATION LOOP ------------------- C----------------------------------------------------------------------- C----------------------------------------------------------------------- Cmp FIM=0. Cmp DO 400 L=LM,1,-1 C----------------------------------------------------------------------- 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,MYPE,L write(6,*) 'Q,T: ', Q(I,J,L),T(I,J,L),MYPE,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 C----------------------------------------------------------------------- !$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 C !$omp parallel do 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 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 C--------------DIAGONAL CONTRIBUTIONS TO PRESSURE GRADIENT FORCE-------- C CC 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 C !$omp parallel do 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) PSE(I,J)=(FIM (I+IHE(J),J-1)-FIM (I,J)) 1 *(DWDT(I+IHE(J),J-1,L)+DWDT(I,J,L)) PPSE(I,J)=PSE(I,J)*ADPDSE(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 C--------------LAT & LONG PRESSURE FORCE COMPONENTS--------------------- C !$omp parallel do private (dcnek,dcsek,dpnek,dpsek) 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 C--------------UPDATE U AND V (CORIOLIS & PGF)-------------------------- C !$omp parallel do private (dpfnek,dpfsek) 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.)-UTK) 1 *VM(I,J)+UTK V(I,J,L)=(VPK-F0K*U(I,J,L)-VTK) 1 *VM(I,J)+VTK 300 CONTINUE C----------------------------------------------------------------------- 400 CONTINUE C----------------------------------------------------------------------- RETURN END