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----------------------------------------------------------------------- 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,*) 'same pintlg at different levels: ', MYPE,I,J,L C write(6,*) 'PINT vals: ', PINT(I,J,L),PINT(I,J,L-1), + PINT(I,J,L-2),PINT(I,J,L-3) C 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) CHOU if equal PINTs replace by 4 point average write(6,*) '>>PINT 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) PINT(I,J,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))*0.25 PINTLG(I,J,L-1)=ALOG(PINT(I,J,L-1)) write(6,*) '>>mype,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) 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 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 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