A ROTINA DA ESQUERDA É A VERSÃO ANTIGA DE 2012 E A ROTINA DA DIREITA É A VERSÃO MODULAR DO Eta1km !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& | SUBROUTINE EPS SUBROUTINE EPS | !>---------------------------------------------------------------------------------------------- C ****************************************************************** | !> SUBROUTINE EPS C$$$ SUBPROGRAM DOCUMENTATION BLOCK | !> C . . . | !> SUBPROGRAM: EPS - ????? C SUBPROGRAM: EPS | !> PROGRAMMER: JANJIC C PRGRMMR: JANJIC ORG: W/NP22 DATE: 9?-??-?? | !> ORG: W/NP22 C | !> DATE: 9?-??-?? C ABSTRACT: | !> C EPS COMPUTES THE VERTICAL AND HORIZONTAL ADVECTION OF DZ/DT | !> ABSTRACT: C | !> EPS COMPUTES THE VERTICAL AND HORIZONTAL ADVECTION OF DZ / DT. C PROGRAM HISTORY LOG: | !> C 9?-??-?? JANJIC - ORIGINATOR | !> PROGRAM HISTORY LOG: C 00-01-05 BLACK - DISTRIBUTED MEMORY AND THREADS | !> 9?-??-?? JANJIC - ORIGINATOR C 26-11-12 RISTIC IVAN - LOT OF CHANGES AND BUG FIX | !> 00-01-05 BLACK - DISTRIBUTED MEMORY AND THREADS C | !> 26-11-12 RISTIC IVAN - LOT OF CHANGES AND BUG FIX C USAGE: CALL EPS FROM MAIN PROGRAM | !> 18-01-15 LUCCI - MODERNIZATION OF THE CODE, INCLUDING: C INPUT ARGUMENT LIST: | !> * F77 TO F90/F95 C NONE | !> * INDENTATION & UNIFORMIZATION CODE C | !> * REPLACEMENT OF COMMONS BLOCK FOR MODULES C OUTPUT ARGUMENT LIST: | !> * DOCUMENTATION WITH DOXYGEN C NONE | !> * OPENMP FUNCTIONALITY C | !> C OUTPUT FILES: | !> INPUT ARGUMENT LIST: C NONE | !> NONE C | !> C SUBPROGRAMS CALLED: | !> OUTPUT ARGUMENT LIST: C | !> NONE C UNIQUE: NONE | !> C | !> INPUT/OUTPUT ARGUMENT LIST: C LIBRARY: NONE | !> NONE C | !> C COMMON BLOCKS: CTLBLK | !> USE MODULES: CONTIN C LOOPS | !> CTLBLK C MASKS | !> DYNAM C DYNAM | !> F77KINDS C CONTIN | !> GLB_TABLE C VRBLS | !> INDX C PVRBLS | !> LOOPS C NHYDRO | !> MAPPINGS C INDX | !> MASKS C ATTRIBUTES: | !> MPPCOM C LANGUAGE: FORTRAN 90 | !> NHYDRO C MACHINE : IBM SP | !> PARMETA C$$$ | !> PVRBLS C*********************************************************************** | !> TEMPCOM INCLUDE "EXCHM.h" | !> TOPO INCLUDE "parmeta" | !> VRBLS INCLUDE "mpp.h" | !> INCLUDE "mpif.h" | !> DRIVER : EBU C----------------------------------------------------------------------- | !> P A R A M E T E R | !> CALLS : EXCH &(IM1=IM-1,JAM=6+2*(JM-10) | !> MPI_ALLREDUCE &,JAMD=(JAM*2-10)*3,LP1=LM+1) | !----------------------------------------------------------------------------------------------- C----------------------------------------------------------------------- | USE CONTIN PARAMETER (EPSQ=1.E-12,EPSQ2=0.2) | USE CTLBLK PARAMETER (FF1=0.52500,FF2=-0.64813,FF3=0.24520,FF4=-0.12189) | USE DYNAM C----------------------------------------------------------------------- | USE EXCHM P A R A M E T E R | USE F77KINDS &(NTSHY=3,EPSFC=1./1.05,EPSN=-EPSFC,EPSP=EPSFC,ZERO=1.E-06 | USE GLB_TABLE &,G=9.8,CP=1004.6,CAPA=287.04/CP,GMA=-287.04*(1.-CAPA)/2.) | USE INDX C----------------------------------------------------------------------- | USE LOOPS L O G I C A L | USE MAPOT & RUN,FIRST,RESTRT,SIGMA,TOP,BOT | USE MAPPINGS C----------------------------------------------------------------------- | USE MASKS INCLUDE "CTLBLK.comm" | USE MPPCOM C----------------------------------------------------------------------- | USE NHYDRO INCLUDE "LOOPS.comm" | USE PARMETA C----------------------------------------------------------------------- | USE PVRBLS INCLUDE "MASKS.comm" | USE TEMPCOM C----------------------------------------------------------------------- | USE TOPO INCLUDE "DYNAM.comm" | USE VRBLS C----------------------------------------------------------------------- | ! INCLUDE "CONTIN.comm" | IMPLICIT NONE C----------------------------------------------------------------------- | ! INCLUDE "VRBLS.comm" | INCLUDE "EXCHM.h" C----------------------------------------------------------------------- | INCLUDE "mpif.h" INCLUDE "PVRBLS.comm" | ! C----------------------------------------------------------------------- | INTEGER(KIND=I4KIND), PARAMETER :: IM1 = IM - 1 INCLUDE "NHYDRO.comm" | INTEGER(KIND=I4KIND), PARAMETER :: JAMD = (JAM * 2 - 10) * 3 C----------------------------------------------------------------------- | INTEGER(KIND=I4KIND), PARAMETER :: NTSHY = 3 INCLUDE "INDX.comm" | ! C----------------------------------------------------------------------- | REAL (KIND=R4KIND), PARAMETER :: EPSQ = 1.E-12 C-------------LOCAL VARIABLES REQUIRING DECLARATION--------------------- | REAL (KIND=R4KIND), PARAMETER :: EPSQ2 = 0.2 C----------------------------------------------------------------------- | REAL (KIND=R4KIND), PARAMETER :: FF1 = 0.52500 R E A L | REAL (KIND=R4KIND), PARAMETER :: FF2 = -0.64813 & PONE(LM+1),PSTR(LM+1),PNP1(LM+1),COFF(LM+1),DFRC(LM+1) | REAL (KIND=R4KIND), PARAMETER :: FF3 = 0.24520 &,CHI(LM+1),CHIN(LM+1),WRES(LM+1),CLIM(LM+1) | REAL (KIND=R4KIND), PARAMETER :: FF4 = -0.12189 R E A L | REAL (KIND=R4KIND), PARAMETER :: EPSFC = 1. / 1.05 & RDPP(LM),C0(LM),B1(LM),B2(LM),B3(LM) | REAL (KIND=R4KIND), PARAMETER :: EPSN = -EPSFC C----------------------------------------------------------------------- | REAL (KIND=R4KIND), PARAMETER :: EPSP = EPSFC REAL W3(LM),W4(LM) | REAL (KIND=R4KIND), PARAMETER :: ZERO = 1.E-06 REAL ETADTL(LM),DWL(LM),AFR(LM) | REAL (KIND=R4KIND), PARAMETER :: G = 9.8 INTEGER LA(LM) | REAL (KIND=R4KIND), PARAMETER :: CP = 1004.6 C----------------------------------------------------------------------- | REAL (KIND=R4KIND), PARAMETER :: CAPA = 287.04 / CP R E A L | REAL (KIND=R4KIND), PARAMETER :: GMA = -287.04 * (1. - CAPA) / 2. & AFP (IDIM1:IDIM2,JDIM1:JDIM2,LM),AFQ (IDIM1:IDIM2,JDIM1:JDIM2,LM) | ! &,W1 (IDIM1:IDIM2,JDIM1:JDIM2,LM),DWST(IDIM1:IDIM2,JDIM1:JDIM2,LM) | LOGICAL(KIND=L4KIND) &,DARE(IDIM1:IDIM2,JDIM1:JDIM2) ,DVOL(IDIM1:IDIM2,JDIM1:JDIM2,LM) | & TOP , BOT &,EMH (IDIM1:IDIM2,JDIM1:JDIM2) | !------------------------ &,ANE (IDIM1:IDIM2,JDIM1:JDIM2) ,ASE (IDIM1:IDIM2,JDIM1:JDIM2) | ! IMPLICIT NONE VARIABLES C----------------------------------------------------------------------- | !------------------------ I N T E G E R | INTEGER(KIND=I4KIND) & IFPA(IDIM1:IDIM2,JDIM1:JDIM2,LM),IFQA(IDIM1:IDIM2,JDIM1:JDIM2,LM) | & I , J , K , LMP , LAP , LLAP , JFP , JFQ , IRECV , &,IFPF(IDIM1:IDIM2,JDIM1:JDIM2,LM),IFQF(IDIM1:IDIM2,JDIM1:JDIM2,LM) | & JHL , JHH , IHL , IHH , JX , IX , KN , KP &,JFPA(IDIM1:IDIM2,JDIM1:JDIM2,LM),JFQA(IDIM1:IDIM2,JDIM1:JDIM2,LM) | ! &,JFPF(IDIM1:IDIM2,JDIM1:JDIM2,LM),JFQF(IDIM1:IDIM2,JDIM1:JDIM2,LM) | REAL (KIND=R4KIND) C----------------------------------------------------------------------- | & ADDT , RDT , RR , ARR , DWP , DETAL , RFACW , W4P , AFRP , R E A L | & D2PQW , WP , W00 , WP0 , ENH , HM , DVOLP , TTA , TTB , & GSUMS(2,LM),XSUMS(2,LM),SUMPW,SUMNW | & PP , QP , DWSTIJ , RFWIJ , W1IJ , WSTIJ , W0Q , WA , DWDTMX , C----------------------------------------------------------------------- | & DWDTMN , DWDTT , GDT , GDT2 , WGHT , FFC , PDP , RDP , DPPL , C*********************************************************************** | & DPSTR , PP1 , TFC , TTFC , PSTRUP , RDPDN , RDPUP , RPD , PSTRDN , C----------------------------------------------------------------------- | & TMP , HBM2IJ , DPTU , FCT , DPTL , DELP IF(NTSD.LE.NTSHY.OR.HYDRO)THEN | !-------------------------------------- C*** | ! LOCAL VARIABLES REQUIRING DECLARATION !$omp parallel do | !-------------------------------------- DO J=MYJS_P2,MYJE_P2 | REAL (KIND=R4KIND), DIMENSION(LM+1) DO I=MYIS_P1,MYIE_P1 | & PONE , PSTR , PNP1 , COFF , DFRC , CHI , CHIN , WRES , CLIM PINT(I,J,1)=PT | ! ENDDO | REAL (KIND=R4KIND), DIMENSION(LM) ENDDO | & RDPP , C0 , B1 , B2 , B3 , W3 , W4 C | ! DO L=1,LM | REAL (KIND=R4KIND), DIMENSION(LM) !$omp parallel do | & ETADTL , DWL , AFR DO J=MYJS_P2,MYJE_P2 | ! DO I=MYIS_P1,MYIE_P1 | INTEGER(KIND=I4KIND), DIMENSION(LM) DWDT(I,J,L)=1. | & LA PDWDT(I,J,L)=1. | ! PINT(I,J,L+1)=PDSL(I,J)*DETA(L)+PINT(I,J,L) | REAL (KIND=R4KIND), DIMENSION(IDIM1:IDIM2, JDIM1:JDIM2, LM) ENDDO | & AFP , AFQ , W1 , DWST , DVOL ENDDO | ! ENDDO | REAL (KIND=R4KIND), DIMENSION(IDIM1:IDIM2, JDIM1:JDIM2) C*** | & DARE , EMH , ANE , ASE RETURN | ! C*** | INTEGER(KIND=I4KIND), DIMENSION(IDIM1:IDIM2, JDIM1:JDIM2, LM) ENDIF | & IFPA , IFQA , C----------------------------------------------------------------------- | & IFPF , IFQF , ADDT=DT | & JFPA , JFQA , RDT=1./ADDT | & JFPF , JFQF C-----------------TIME TENDENCY----------------------------------------- | ! !$omp parallel do | REAL (KIND=R4KIND), DIMENSION(2, LM) DO L=1,LM | & GSUMS , XSUMS DO J=MYJS_P1,MYJE_P1 | ! DO I=MYIS_P1,MYIE_P1 | REAL (KIND=R4KIND) DWDT(I,J,L)=(W(I,J,L)-DWDT(I,J,L))*HTM(I,J,L)*HBM2(I,J)*RDT | & SUMPW , SUMNW ENDDO | ! ENDDO | IF (NTSD <= NTSHY .OR. HYDRO) THEN ENDDO | !------- C | ! OPENMP C----------------------------------------------------------------------- | !------- C-----------------VERTICAL ADVECTION------------------------------------ | ! C----------------------------------------------------------------------- | !$omp parallel do C | ! !$omp parallel do | DO J=MYJS_P2,MYJE_P2 !$omp& private (afr,afrp,arr,bot,d2pqw,detal,dwl,dwp,etadtl, | DO I=MYIS_P1,MYIE_P1 !$omp& la,lap,llap,lmp,rfacw,rr,sumnw,sumpw,top, | PINT(I,J,1) = PT !$omp& w00,w3,w4,w4p,wp,wp0) | END DO DO 200 J=MYJS_P1,MYJE_P1 | END DO DO 200 I=MYIS_P1,MYIE_P1 | C----------------------------------------------------------------------- | DO K=1,LM LMP=LMH(I,J) | !------- C----------------------------------------------------------------------- | ! OPENMP DO L=1,LMP | !------- W3(L)=W(I,J,L) | ! W4(L)=W3(L) | !$omp parallel do ENDDO | ! C | DO J=MYJS_P2,MYJE_P2 ETADTL(1)=ETADT(I,J,1)*0.5 | DO I=MYIS_P1,MYIE_P1 C | DWDT(I,J,K) = 1. DO L=2,LMP-1 | PDWDT(I,J,K) = 1. ETADTL(L)=(ETADT(I,J,L-1)+ETADT(I,J,L))*0.5 | PINT(I,J,K+1) = PDSL(I,J) * DETA(K) + PINT(I,J,K) ENDDO | END DO C | ENDDO ETADTL(LMP)=ETADT(I,J,LMP-1)*0.5 | END DO C----------------------------------------------------------------------- | ! SUMPW=0. | RETURN SUMNW=0. | ! C | END IF DO L=1,LMP | ! RR=ETADTL(L)*(-ADDT) | ADDT = DT ARR=ABS(RR) | RDT = 1. / ADDT C | !-------------- IF(ARR.GT.0.)THEN | ! TIME TENDENCY LAP=RR/ARR | !-------------- ELSE | !------- LAP=0 | ! OPENMP ENDIF | !------- C | ! LA(L)=LAP | !$omp parallel do LLAP=L+LAP | ! C | DO K=1,LM TOP=.FALSE. | DO J=MYJS_P1,MYJE_P1 BOT=.FALSE. | DO I=MYIS_P1,MYIE_P1 C | DWDT(I,J,K) = (W(I,J,K) - DWDT(I,J,K)) * HTM(I,J,K) * HBM2(I,J) * RDT IF(LLAP.GT.0.AND.LLAP.LT.LMH(I,J)+1.AND.LAP.NE.0)THEN | END DO RR=ARR/ABS(AETA(LLAP)-AETA(L)) | END DO AFR(L)=(((FF4*RR+FF3)*RR+FF2)*RR+FF1)*RR | END DO DWP=(W3(LLAP)-W3(L))*RR | !------------------- DWL(L)=DWP | ! VERTICAL ADVECTION ELSE | !------------------- TOP=LLAP.EQ.0 | !------- BOT=LLAP.EQ.LMH(I,J)+1 | ! OPENMP C | !------- RR=0. | ! AFR(L)=0. | !$omp parallel do DWL(L)=0. | !$omp private (AFR , AFRP , ARR , BOT , D2PQW , DETAL , DWL , DWP , ENDIF | !$omp ETADTL , LA , LAP , LLAP , LMP , RFACW , RR , SUMNW , ENDDO | !$omp SUMPW , TOP , W00 , W3 , W4 , W4P , WP , WP0 ) C | ! IF(TOP)THEN | DO 200 J=MYJS_P1,MYJE_P1 IF(LA(2).LT.0)THEN | DO 200 I=MYIS_P1,MYIE_P1 DWL(1)=-DWL(2)*DETA(2)/DETA(1) | ! ENDIF | LMP=LMH(I,J) ENDIF | ! C | DO K=1,LMP IF(BOT)THEN | W3(K) = W(I,J,K) IF(LA(LMP-1).GT.0)THEN | W4(K) = W3(K) DWL(LMP)=-DWL(LMP-1)*DETA(LMP-1)/DETA(LMP) | END DO ENDIF | ENDIF | ETADTL(1) = ETADT(I,J,1) * 0.5 C | DO L=1,LMP | DO K=2,LMP-1 DETAL=DETA(L) | ETADTL(K) = (ETADT(I,J,K-1) + ETADT(I,J,K)) * 0.5 DWP=DWL(L)*DETAL | END DO IF(DWP.GT.0.)THEN | ! SUMPW=SUMPW+DWP | ETADTL(LMP) = ETADT(I,J,LMP-1) * 0.5 ELSE | ! SUMNW=SUMNW+DWP | SUMPW = 0. ENDIF | SUMNW = 0. ENDDO | C--------------FIRST MOMENT CONSERVING FACTOR--------------------------- | DO K=1,LMP IF(SUMPW.GT.1.E-9)THEN | RR = ETADTL(K) * (-ADDT) RFACW=-SUMNW/SUMPW | ARR = ABS(RR) ELSE | RFACW=1. | IF (ARR > 0.) THEN ENDIf | LAP = RR / ARR C | ELSE IF(RFACW.LT.0.9.OR.RFACW.GT.1.1)RFACW=1. | LAP = 0 C--------------IMPOSE CONSERVATION ON ADVECTION------------------------- | END IF IF(RFACW.LT.1.)THEN | ! DO L=1,LMP | LA(K) = LAP DWP=DWL(L) | LLAP = K + LAP IF(DWP.LT.0.)DWP=DWP/RFACW | ! W4(L)=W3(L)+DWP | TOP = .FALSE. ENDDO | BOT = .FALSE. ELSE | ! DO L=1,LMP | IF (LLAP > 0 .AND. LLAP < LMH(I,J)+1 .AND. LAP /= 0) THEN DWP=DWL(L) | RR = ARR / ABS(AETA(LLAP) - AETA(K)) IF(DWP.GE.0.)DWP=DWP*RFACW | AFR(K) = (((FF4 * RR + FF3) * RR + FF2) * RR + FF1) * RR W4(L)=W3(L)+DWP | DWP = (W3(LLAP) - W3(K)) * RR ENDDO | DWL(K) = DWP ENDIF | ELSE C--------------ANTI-FILTERING STEP-------------------------------------- | TOP = LLAP == 0 SUMPW=0. | BOT = LLAP == LMH(I,J)+1 SUMNW=0. | ! C--------------ANTI-FILTERING LIMITERS---------------------------------- | RR = 0. DO 50 L=2,LMP-1 | AFR(K) = 0. DETAL=DETA(L) | DWL(K) = 0. C | END IF W4P=W4(L) | END DO C | ! LAP=LA(L) | IF (TOP) THEN C | IF (LA(2) < 0) THEN IF(LAP.NE.0)THEN | DWL(1) = -DWL(2) * DETA(2) / DETA(1) AFRP=2.*AFR(L)*(AETA(L+LAP)-AETA(L))**2 | END IF 1 /(AETA(L+LAP)-AETA(L-LAP)) | END IF D2PQW=((W4(L+LAP)-W4P)/(AETA(L+LAP)-AETA(L)) | ! 1 -(W4P-W4(L-LAP))/(AETA(L)-AETA(L-LAP)))*AFRP | IF (BOT) THEN ELSE | IF (LA(LMP-1) > 0)THEN D2PQW=0. | DWL(LMP) = -DWL(LMP-1) * DETA(LMP-1) / DETA(LMP) ENDIF | END IF C | END IF WP=W4P-D2PQW | ! C | DO K=1,LMP W00=W3(L) | DETAL = DETA(K) WP0=W3(L+LAP) | DWP = DWL (K) * DETAL C | IF (DWP > 0.) THEN WP=MAX(WP,MIN(W00,WP0)) | SUMPW = SUMPW + DWP WP=MIN(WP,MAX(W00,WP0)) | ELSE C | SUMNW = SUMNW + DWP DWP=WP-W4P | END IF C | END DO DWL(L)=DWP | !------------------------------- C | ! FIRST MOMENT CONSERVING FACTOR DWP=DWP*DETAL | !------------------------------- C | IF (SUMPW > 1.E-9) THEN IF(DWP.GT.0.)THEN | RFACW = -SUMNW / SUMPW SUMPW=SUMPW+DWP | ELSE ELSE | RFACW = 1. SUMNW=SUMNW+DWP | END IF ENDIF | ! C | IF(RFACW < 0.9 .OR. RFACW > 1.1) RFACW = 1. 50 CONTINUE | !--------------------------------- C----------------------------------------------------------------------- | ! IMPOSE CONSERVATION ON ADVECTION DWL(1)=0. | !--------------------------------- C | IF (RFACW < 1.) THEN DWL(LMP)=0. | DO K=1,LMP C--------------FIRST MOMENT CONSERVING FACTOR--------------------------- | DWP = DWL(K) IF(SUMPW.GT.1.E-9)THEN | IF (DWP < 0.) DWP = DWP / RFACW RFACW=-SUMNW/SUMPW | W4(K) = W3(K) + DWP ELSE | END DO RFACW=1. | ELSE ENDIF | DO K=1,LMP C | DWP = DWL(K) IF(RFACW.LT.0.9.OR.RFACW.GT.1.1)RFACW=1. | IF (DWP >= 0.) DWP = DWP * RFACW C--------------IMPOSE CONSERVATION ON ANTI-FILTERING-------------------- | W4(K) = W3(K) + DWP IF(RFACW.LT.1.)THEN | END DO DO L=1,LMP | END IF DWP=DWL(L) | !-------------------- IF(DWP.GE.0.)DWP=DWP*RFACW | ! ANTI-FILTERING STEP DWDT(I,J,L)=DWDT(I,J,L)-(W4(L)+DWP-W3(L))*RDT | !-------------------- ENDDO | SUMPW = 0. ELSE | SUMNW = 0. DO L=1,LMP | !------------------------ DWP=DWL(L) | ! ANTI-FILTERING LIMITERS IF(DWP.LT.0.)DWP=DWP/RFACW | !------------------------ DWDT(I,J,L)=DWDT(I,J,L)-(W4(L)+DWP-W3(L))*RDT | DO 50 K=2,LMP-1 ENDDO | DETAL = DETA(K) ENDIF | ! C----------------------------------------------------------------------- | W4P = W4(K) 200 CONTINUE | ! C--------------END OF VERTICAL ADVECTION-------------------------------- | LAP = LA(K) C | ! C----------------------------------------------------------------------- | IF (LAP /= 0) THEN C*********************************************************************** | AFRP = 2. * AFR(K) * (AETA(K+LAP) - AETA(K)) ** 2 / (AETA(K+LAP) - AETA(K-L ENH=ADDT/(08.*DY) | D2PQW = ((W4(K+LAP) - W4P) / (AETA(K+LAP) - AETA(K)) C----------------------------------------------------------------------- | & - (W4P - W4(K-LAP)) / (AETA(K) - AETA(K-LAP))) !$omp parallel do | & * AFRP DO J=MYJS_P2,MYJE_P2 | ELSE DO I=MYIS_P1,MYIE_P1 | D2PQW = 0. EMH(I,J)=ADDT/(08.*DX(I,J)) | END IF DARE(I,J)=HBM2(I,J)*DX(I,J)*DY | ! ENDDO | WP = W4P - D2PQW ENDDO | ! C----------------------------------------------------------------------- | W00 = W3(K) C---------------------HORIZONTAL ADVECTION------------------------------ | WP0 = W3(K+LAP) C----------------------------------------------------------------------- | ! C*********************************************************************** | WP = MAX(WP,MIN(W00,WP0)) !$omp parallel do | WP = MIN(WP,MAX(W00,WP0)) !$omp& private (dvolp,dwstij,hm,jfp,jfq,pp,qp, | ! !$omp& sumnw,sumpw,tta,ttb) | DWP = WP - W4P DO 300 L=1,LM | ! C*********************************************************************** | DWL(K) = DWP DO J=MYJS_P2,MYJE_P2 | ! DO I=MYIS_P1,MYIE_P1 | DWP = DWP * DETAL DVOL(I,J,L)=DARE(I,J)*PDSL(I,J)*DETA(L) | ! HM=HTM(I,J,L) | IF (DWP > 0.) THEN W1(I,J,L)=W(I,J,L)*HM | SUMPW = SUMPW + DWP ENDDO | ELSE ENDDO | SUMNW = SUMNW + DWP C----------------------------------------------------------------------- | END IF SUMPW=0. | ! SUMNW=0. | 50 END DO C | ! DO 225 J=MYJS2_P1,MYJE2_P1 | DWL(1) = 0. DO 225 I=MYIS1_P1,MYIE1_P1 | C | DWL(LMP) = 0. DVOLP=DVOL(I,J,L)*HBM3(I,J) | !------------------------------- TTA=(U(I,J-1,L)+U(I+IHW(J),J,L)+U(I+IHE(J),J,L)+U(I,J+1,L)) | ! FIRST MOMENT CONSERVING FACTOR 2 *HBM2(I,J)*EMH(I,J) | !------------------------------- TTB=(V(I,J-1,L)+V(I+IHW(J),J,L)+V(I+IHE(J),J,L)+V(I,J+1,L)) | IF (SUMPW > 1.E-9) THEN 2 *HBM2(I,J)*ENH | RFACW = -SUMNW / SUMPW C | ELSE PP=-TTA-TTB | RFACW = 1. QP= TTA-TTB | END IF C | ! JFP=INT(SIGN(1.,PP)) | IF (RFACW < 0.9 .OR. RFACW > 1.1) RFACW = 1. JFQ=INT(SIGN(1.,QP)) | !-------------------------------------- C | ! IMPOSE CONSERVATION ON ANTI-FILTERING IFPA(I,J,L)=IHE(J)+I+( JFP-1)/2 | !-------------------------------------- IFQA(I,J,L)=IHE(J)+I+(-JFQ-1)/2 | IF (RFACW < 1.) THEN C | DO K=1,LMP JFPA(I,J,L)=J+JFP | DWP = DWL(K) JFQA(I,J,L)=J+JFQ | IF (DWP >= 0.) DWP = DWP * RFACW C | DWDT(I,J,K) = DWDT(I,J,K) - (W4(K) + DWP - W3(K)) * RDT IFPF(I,J,L)=IHE(J)+I+(-JFP-1)/2 | END DO IFQF(I,J,L)=IHE(J)+I+( JFQ-1)/2 | ELSE C | DO K=1,LMP JFPF(I,J,L)=J-JFP | DWP = DWL(K) JFQF(I,J,L)=J-JFQ | IF (DWP < 0.) DWP = DWP / RFACW C | DWDT(I,J,K) = DWDT(I,J,K) - (W4(K) + DWP - W3(K)) * RDT PP=ABS(PP)*HTM(I,J,L)*HTM(IFPA(I,J,L),JFPA(I,J,L),L) | END DO QP=ABS(QP)*HTM(I,J,L)*HTM(IFQA(I,J,L),JFQA(I,J,L),L) | END IF C | ! AFP(I,J,L)=(((FF4*PP+FF3)*PP+FF2)*PP+FF1)*PP | 200 END DO AFQ(I,J,L)=(((FF4*QP+FF3)*QP+FF2)*QP+FF1)*QP | !-------------------------- C | ! END OF VERTICAL ADVECTION DWSTIJ=(W(IFPA(I,J,L),JFPA(I,J,L),L)-W(I,J,L))*PP | !-------------------------- 2 +(W(IFQA(I,J,L),JFQA(I,J,L),L)-W(I,J,L))*QP | ! C | ENH = ADDT / (08. * DY) DWST(I,J,L)=DWSTIJ | !------- C | ! OPENMP 225 CONTINUE | !------- C*** | ! C*** GLOBAL SUM FOR CONSERVATION | !$omp parallel do C*** | ! DO 230 J=MYJS2,MYJE2 | DO J=MYJS_P2,MYJE_P2 DO 230 I=MYIS1,MYIE1 | DO I=MYIS_P1,MYIE_P1 C | EMH(I,J) = ADDT / (08. * DX(I,J)) DVOLP=DVOL(I,J,L)*HBM3(I,J) | DARE(I,J) = HBM2(I,J) * DX(I,J) * DY DWSTIJ=DWST(I,J,L)*DVOLP | END DO C | END DO IF(DWSTIJ.GT.0.)THEN | !--------------------- SUMPW=SUMPW+DWSTIJ | ! HORIZONTAL ADVECTION ELSE | !--------------------- SUMNW=SUMNW+DWSTIJ | !------- ENDIF | ! OPENMP C | !------- 230 CONTINUE | ! C | !$omp parallel do C----------------------------------------------------------------------- | !$omp private (DVOLP , DWSTIJ , HM , JFP , JFQ , PP , QP , SUMNW , XSUMS(1,L)=SUMPW | !$omp SUMPW , TTA , TTB ) XSUMS(2,L)=SUMNW | ! C | DO 300 K=1,LM 300 CONTINUE ! End of Vertical Loop | ! C----------------------------------------------------------------------- | DO J=MYJS_P2,MYJE_P2 C | DO I=MYIS_P1,MYIE_P1 C*** GLOBAL REDUCTION | DVOL(I,J,K) = DARE(I,J) * PDSL(I,J) * DETA(K) C | HM = HTM(I,J,K) CALL MPI_ALLREDUCE(XSUMS,GSUMS,2*LM,MPI_REAL,MPI_SUM, | W1(I,J,K) = W(I,J,K) * HM 1 MPI_COMM_COMP,IRECV) | END DO IF(IRECV.NE.0)THEN | END DO PRINT*, ' RETURN CODE FROM 1st ALLREDUCE IN EPS = ',IRECV | ! STOP | SUMPW = 0. ENDIF | SUMNW = 0. C----------------------------------------------------------------------- | ! !$omp parallel do | DO 225 J=MYJS2_P1,MYJE2_P1 !$omp& private(d2pqw,dvolp,dwstij,rfacw,rfwij,sumnw,sumpw, | DO 225 I=MYIS1_P1,MYIE1_P1 !$omp& w00,w0q,w1ij,wp0,wstij) | ! C----------------------------------------------------------------------- | DVOLP = DVOL(I,J,K) * HBM3(I,J) DO 400 L=1,LM | TTA = (U(I,J-1,K) + U(I+IHW(J),J,K) + U(I+IHE(J),J,K) + U(I,J+1,K)) C----------------------------------------------------------------------- | & * HBM2(I,J) * EMH(I,J) C | ! SUMPW=GSUMS(1,L) | TTB = (V(I,J-1,K) + V(I+IHW(J),J,K) + V(I+IHE(J),J,K) + V(I,J+1,K)) SUMNW=GSUMS(2,L) | & * HBM2(I,J) * ENH C | ! C--------------FIRST MOMENT CONSERVING FACTOR--------------------------- | PP = -TTA - TTB IF(SUMPW.GT.1.)THEN | QP = TTA - TTB RFACW=-SUMNW/SUMPW | ! ELSE | JFP = INT(SIGN(1.,PP)) RFACW=1. | JFQ = INT(SIGN(1.,QP)) ENDIF | ! IF(RFACW.LT.0.9.OR.RFACW.GT.1.1)RFACW=1. | IFPA(I,J,K) = IHE(J) + I + ( JFP-1) / 2 C--------------IMPOSE CONSERVATION ON ADVECTION------------------------- | IFQA(I,J,K) = IHE(J) + I + (-JFQ-1) / 2 IF(RFACW.LT.1.)THEN | ! DO J=MYJS2_P1,MYJE2_P1 | JFPA(I,J,K) = J + JFP DO I=MYIS1_P1,MYIE1_P1 | JFQA(I,J,K) = J + JFQ DWSTIJ=DWST(I,J,L) | ! RFWIJ=HBM3(I,J)*(RFACW-1.)+1. | IFPF(I,J,K) = IHE(J) + I + (-JFP-1) / 2 IF(DWSTIJ.LT.0.)DWSTIJ=DWSTIJ/RFWIJ | IFQF(I,J,K) = IHE(J) + I + ( JFQ-1) / 2 W1(I,J,L)=W(I,J,L)+DWSTIJ | ! ENDDO | JFPF(I,J,K) = J - JFP ENDDO | JFQF(I,J,K) = J - JFQ ELSE | ! DO J=MYJS2_P1,MYJE2_P1 | PP = ABS(PP) * HTM(I,J,K) * HTM(IFPA(I,J,K), JFPA(I,J,K),K) DO I=MYIS1_P1,MYIE1_P1 | QP = ABS(QP) * HTM(I,J,K) * HTM(IFQA(I,J,K), JFQA(I,J,K),K) DWSTIJ=DWST(I,J,L) | ! RFWIJ=HBM3(I,J)*(RFACW-1.)+1. | AFP(I,J,K) = (((FF4 * PP + FF3) * PP + FF2) * PP + FF1) * PP IF(DWSTIJ.GE.0.)DWSTIJ=DWSTIJ*RFWIJ | AFQ(I,J,K) = (((FF4 * QP + FF3) * QP + FF2) * QP + FF1) * QP W1(I,J,L)=W(I,J,L)+DWSTIJ | ! ENDDO | DWSTIJ = (W(IFPA(I,J,K), JFPA(I,J,K),K) - W(I,J,K)) * PP ENDDO | & + (W(IFQA(I,J,K), JFQA(I,J,K),K) - W(I,J,K)) * QP ENDIF | ! C--------------ANTI-FILTERING STEP-------------------------------------- | DWST(I,J,K) = DWSTIJ SUMPW=0. | ! SUMNW=0. | 225 END DO C--------------ANTI-FILTERING LIMITERS---------------------------------- | !---------------------------- C | ! GLOBAL SUM FOR CONSERVATION DO 350 J=MYJS2,MYJE2 | !---------------------------- DO 350 I=MYIS1,MYIE1 | DO 230 J=MYJS2,MYJE2 C | DO 230 I=MYIS1,MYIE1 DVOLP=DVOL(I,J,L) | ! W1IJ =W1(I,J,L) | DVOLP = DVOL(I,J,K) * HBM3(I,J) C | DWSTIJ = DWST(I,J,K) * DVOLP D2PQW=((W1(IFPA(I,J,L),JFPA(I,J,L),L)-W1IJ ) | ! 2 -(W1IJ -W1(IFPF(I,J,L),JFPF(I,J,L),L)) | IF (DWSTIJ > 0.) THEN 3 *HTM(IFPF(I,J,L),JFPF(I,J,L),L))*AFP(I,J,L) | SUMPW = SUMPW + DWSTIJ 4 +((W1(IFQA(I,J,L),JFQA(I,J,L),L)-W1IJ ) | ELSE 5 -(W1IJ -W1(IFQF(I,J,L),JFQF(I,J,L),L)) | SUMNW = SUMNW + DWSTIJ 6 *HTM(IFQF(I,J,L),JFQF(I,J,L),L))*AFQ(I,J,L) | END IF C | ! WSTIJ=W1IJ-D2PQW | 230 END DO C | ! W00=W(I ,J ,L) | XSUMS(1,K) = SUMPW WP0=W(IFPA(I,J,L),JFPA(I,J,L),L) | XSUMS(2,K) = SUMNW W0Q=W(IFQA(I,J,L),JFQA(I,J,L),L) | C | CONTINUE !END OF VERTICAL LOOP WSTIJ=AMAX1(WSTIJ,AMIN1(W00,WP0,W0Q)) | ! WSTIJ=AMIN1(WSTIJ,AMAX1(W00,WP0,W0Q)) | 300 END DO C | !----------------- DWSTIJ=WSTIJ-W1IJ | ! GLOBAL REDUCTION C | !----------------- DWST(I,J,L)=DWSTIJ | CALL MPI_ALLREDUCE(XSUMS, GSUMS, 2*LM, MPI_REAL, MPI_SUM, MPI_COMM_COMP, IRECV) C | ! DWSTIJ=DWSTIJ*DVOLP | IF (IRECV /= 0) THEN C | PRINT*, ' RETURN CODE FROM 1ST ALLREDUCE IN EPS = ', IRECV IF(DWSTIJ.GT.0.)THEN | STOP SUMPW=SUMPW+DWSTIJ | END IF ELSE | !------- SUMNW=SUMNW+DWSTIJ | ! OPENMP ENDIF | !------- C | ! 350 CONTINUE | !$omp parallel do C----------------------------------------------------------------------- | !$omp private (D2PQW , DVOLP , DWSTIJ , RFACW , RFWIJ , SUMNW , SUMPW , W00 , XSUMS(1,L)=SUMPW | !$omp W0Q , W1IJ , WP0 , WSTIJ ) XSUMS(2,L)=SUMNW | ! C | DO 400 K=1,LM 400 CONTINUE ! End of Vertical Loop | ! C----------------------------------------------------------------------- | SUMPW = GSUMS(1,K) C | SUMNW = GSUMS(2,K) C*** GLOBAL REDUCTION | !------------------------------- C | ! FIRST MOMENT CONSERVING FACTOR CALL MPI_ALLREDUCE(XSUMS,GSUMS,2*LM,MPI_REAL,MPI_SUM, | !------------------------------- 1 MPI_COMM_COMP,IRECV) | IF (SUMPW > 1.) THEN IF(IRECV.NE.0)THEN | RFACW = -SUMNW / SUMPW PRINT*, ' RETURN CODE FROM 2nd ALLREDUCE IN EPS = ',IRECV | ELSE STOP | RFACW = 1. ENDIF | END IF C----------------------------------------------------------------------- | ! C | IF (RFACW < 0.9 .OR. RFACW > 1.1) RFACW = 1. !$omp parallel do | !--------------------------------- !$omp& private(dwstij,rfacw,rfwij,sumnw,sumpw) | ! IMPOSE CONSERVATION ON ADVECTION C----------------------------------------------------------------------- | !--------------------------------- DO 425 L=1,LM | IF (RFACW < 1.) THEN C | DO J=MYJS2_P1,MYJE2_P1 SUMPW=GSUMS(1,L) | DO I=MYIS1_P1,MYIE1_P1 SUMNW=GSUMS(2,L) | DWSTIJ = DWST(I,J,K) C | RFWIJ = HBM3(I,J) * (RFACW-1.) + 1. C--------------FIRST MOMENT CONSERVING FACTOR--------------------------- | IF (DWSTIJ < 0.) DWSTIJ = DWSTIJ / RFWIJ IF(SUMPW.GT.1.)THEN | W1(I,J,K) = W(I,J,K) + DWSTIJ RFACW=-SUMNW/SUMPW | END DO ELSE | END DO RFACW=1. | ELSE ENDIF | DO J=MYJS2_P1,MYJE2_P1 IF(RFACW.LT.0.9.OR.RFACW.GT.1.1)RFACW=1. | DO I=MYIS1_P1,MYIE1_P1 C--------------IMPOSE CONSERVATION ON ANTI-FILTERING-------------------- | DWSTIJ = DWST(I,J,K) IF(RFACW.LT.1.)THEN | RFWIJ = HBM3(I,J) * (RFACW-1.) + 1. DO J=MYJS2,MYJE2 | IF (DWSTIJ >= 0.) DWSTIJ = DWSTIJ * RFWIJ DO I=MYIS1,MYIE1 | W1(I,J,K) = W(I,J,K) + DWSTIJ DWSTIJ=DWST(I,J,L) | END DO RFWIJ=HBM2(I,J)*(RFACW-1.)+1. | END DO IF(DWSTIJ.GE.0.)DWSTIJ=DWSTIJ*RFWIJ | END IF DWDT(I,J,L)=(DWDT(I,J,L)-(W1(I,J,L)+DWSTIJ-W(I,J,L))*RDT) | !-------------------- 1 *HBM3(I,J) | ! ANTI-FILTERING STEP ENDDO | !-------------------- ENDDO | SUMPW = 0. ELSE | SUMNW = 0. DO J=MYJS2,MYJE2 | !------------------------ DO I=MYIS1,MYIE1 | ! ANTI-FILTERING LIMITERS DWSTIJ=DWST(I,J,L) | !------------------------ RFWIJ=HBM2(I,J)*(RFACW-1.)+1. | DO 350 J=MYJS2,MYJE2 IF(DWSTIJ.LT.0.)DWSTIJ=DWSTIJ/RFWIJ | DO 350 I=MYIS1,MYIE1 DWDT(I,J,L)=(DWDT(I,J,L)-(W1(I,J,L)+DWSTIJ-W(I,J,L))*RDT) | ! 1 *HBM3(I,J) | DVOLP = DVOL(I,J,K) ENDDO | W1IJ = W1(I,J,K) ENDDO | ! ENDIF | D2PQW = (( W1(IFPA(I,J,K), JFPA(I,J,K), K) - W1IJ) C----------------------------------------------------------------------- | & - (W1IJ - W1(IFPF(I,J,K), JFPF(I,J,K), K)) C*********************************************************************** | & * HTM(IFPF(I,J,K), JFPF(I,J,K), K)) 425 CONTINUE | & * AFP(I,J,K) + ((W1(IFQA(I,J,K), JFQA(I,J,K), K) - W1IJ) C*********************************************************************** | & - (W1IJ - W1(IFQF(I,J,K), JFQF(I,J,K), k)) C | & * HTM(IFQF(I,J,K), JFQF(I,J,K), K)) * AFQ(I,J,K) C--------------RESTRICTING THE ACCELERATION ALONG THE BOUNDARIES-------- | ! C | WSTIJ = W1IJ - D2PQW c JHL=07 | ! JHL=0 | W00 = W(I,J,K) C*** | WP0 = W(IFPA(I,J,K), JFPA(I,J,K), K) IF(JHL.GT.1)THEN | W0Q = W(IFQA(I,J,K), JFQA(I,J,K), K) JHH=JM-JHL+1 | ! C | WSTIJ = AMAX1(WSTIJ, AMIN1(W00,WP0,W0Q)) IHL=JHL/2+1 | WSTIJ = AMIN1(WSTIJ, AMAX1(W00,WP0,W0Q)) IHH=IM-IHL+MOD(J,2) | ! C----------------------------------------------------------------------- | DWSTIJ = WSTIJ - W1IJ !$omp parallel do private (ix,jx) | ! DO 450 L=1,LM | DWST(I,J,K) = DWSTIJ C | ! DO J=1,JHL | DWSTIJ = DWSTIJ * DVOLP IF(J.GE.MY_JS_GLB-JBPAD2.AND.J.LE.MY_JE_GLB+JTPAD2)THEN | ! JX=J-MY_JS_GLB+1 | IF (DWSTIJ > 0.) THEN DO I=1,IM | SUMPW = SUMPW + DWSTIJ IF(I.GE.MY_IS_GLB-ILPAD2.AND.I.LE.MY_IE_GLB+IRPAD2)THEN | ELSE IX=I-MY_IS_GLB+1 | SUMNW = SUMNW + DWSTIJ DWDT(IX,JX,L)=MAX(DWDT(IX,JX,L)*HBM3(IX,JX),-0.03) | END IF DWDT(IX,JX,L)=MIN(DWDT(IX,JX,L)*HBM3(IX,JX), 0.03) | ! ENDIF | 350 END DO ENDDO | ! ENDIF | XSUMS(1,K) = SUMPW ENDDO | XSUMS(2,K) = SUMNW C | ! DO J=JHH,JM | CONTINUE !END OF VERTICAL LOOP IF(J.GE.MY_JS_GLB-JBPAD2.AND.J.LE.MY_JE_GLB+JTPAD2)THEN | 400 END DO JX=J-MY_JS_GLB+1 | !----------------- DO I=1,IM | ! GLOBAL REDUCTION IF(I.GE.MY_IS_GLB-ILPAD2.AND.I.LE.MY_IE_GLB+IRPAD2)THEN | !----------------- IX=I-MY_IS_GLB+1 | CALL MPI_ALLREDUCE(XSUMS, GSUMS, 2*LM, MPI_REAL, MPI_SUM, MPI_COMM_COMP, IRECV) DWDT(IX,JX,L)=MAX(DWDT(IX,JX,L)*HBM3(IX,JX),-0.03) | ! DWDT(IX,JX,L)=MIN(DWDT(IX,JX,L)*HBM3(IX,JX), 0.03) | IF (IRECV /= 0) THEN ENDIF | PRINT*, ' RETURN CODE FROM 2ND ALLREDUCE IN EPS = ', IRECV ENDDO | STOP ENDIF | END IF ENDDO | !------- C | ! OPENMP DO J=1,JM | !------- IF(J.GE.MY_JS_GLB-JBPAD2.AND.J.LE.MY_JE_GLB+JTPAD2)THEN | ! JX=J-MY_JS_GLB+1 | !$omp parallel do DO I=1,IHL | !$omp private (DWSTIJ , RFACW , RFWIJ , SUMNW , SUMPW) IF(I.GE.MY_IS_GLB-ILPAD2.AND.I.LE.MY_IE_GLB+IRPAD2)THEN | ! IX=I-MY_IS_GLB+1 | DO 425 K=1,LM DWDT(IX,JX,L)=MAX(DWDT(IX,JX,L)*HBM3(IX,JX),-0.03) | ! DWDT(IX,JX,L)=MIN(DWDT(IX,JX,L)*HBM3(IX,JX), 0.03) | SUMPW = GSUMS(1,K) ENDIF | SUMNW = GSUMS(2,K) ENDDO | !------------------------------- ENDIF | ! FIRST MOMENT CONSERVING FACTOR ENDDO | !------------------------------- C | IF (SUMPW > 1.) THEN DO J=1,JM | RFACW = -SUMNW / SUMPW IF(J.GE.MY_JS_GLB-JBPAD2.AND.J.LE.MY_JE_GLB+JTPAD2)THEN | ELSE JX=J-MY_JS_GLB+1 | RFACW = 1. DO I=IHH,IM | END IF IF(I.GE.MY_IS_GLB-ILPAD2.AND.I.LE.MY_IE_GLB+IRPAD2)THEN | ! IX=I-MY_IS_GLB+1 | IF (RFACW < 0.9 .OR. RFACW > 1.1) RFACW = 1. DWDT(IX,JX,L)=MAX(DWDT(IX,JX,L)*HBM3(IX,JX),-0.03) | !-------------------------------------- DWDT(IX,JX,L)=MIN(DWDT(IX,JX,L)*HBM3(IX,JX), 0.03) | ! IMPOSE CONSERVATION ON ANTI-FILTERING ENDIF | !-------------------------------------- ENDDO | IF (RFACW < 1.) THEN ENDIF | DO J=MYJS2,MYJE2 ENDDO | DO I=MYIS1,MYIE1 C | DWSTIJ = DWST(I,J,K) 450 CONTINUE | RFWIJ = HBM2(I,J) * (RFACW-1.) + 1. ENDIF | IF (DWSTIJ >= 0.) DWSTIJ = DWSTIJ * RFWIJ C------------------------------------------------------------------------ | DWDT(I,J,K) = (DWDT(I,J,K) - (W1(I,J,K) + DWSTIJ - W(I,J,K)) * RDT) * HBM3(I CALL EXCH(DWDT,LM,1,1) | END DO WA=0.15 | END DO C | ELSE !$omp parallel do private (ane,ase) | DO J=MYJS2,MYJE2 DO 500 L=1,LM | DO I=MYIS1,MYIE1 C | DWSTIJ = DWST(I,J,K) DO J=MYJS1_P1,MYJE2_P1 | RFWIJ = HBM2(I,J) * (RFACW-1.) + 1. DO I=MYIS_P1,MYIE1_P1 | IF (DWSTIJ < 0.) DWSTIJ = DWSTIJ / RFWIJ ANE(I,J)=(DWDT(I+IHE(J),J+1,L)-DWDT(I,J,L)) | DWDT(I,J,K) = (DWDT(I,J,K) - (W1(I,J,K) + DWSTIJ - W(I,J,K)) * RDT) * HBM3(I & *HTM(I,J,L)*HTM(I+IHE(J),J+1,L) | END DO ENDDO | END DO ENDDO | END IF C | ! DO J=MYJS2_P1,MYJE1_P1 | 425 END DO DO I=MYIS_P1,MYIE1_P1 | !-------------------------------------------------- ASE(I,J)=(DWDT(I+IHE(J),J-1,L)-DWDT(I,J,L)) | ! RESTRICTING THE ACCELERATION ALONG THE BOUNDARIES & *HTM(I+IHE(J),J-1,L)*HTM(I,J,L) | !-------------------------------------------------- ENDDO | JHL = 0 ENDDO | ! C | IF (JHL > 1) THEN DO J=MYJS2,MYJE2 | JHH = JM - JHL + 1 DO I=MYIS1,MYIE1 | ! DWDT(I,J,L)=(ANE(I,J)-ANE(I+IHW(J),J-1) | IHL = JHL / 2 + 1 1 +ASE(I,J)-ASE(I+IHW(J),J+1))*WA*HBM2(I,J) | IHH = IM - IHL + MOD(J,2) 2 +DWDT(I,J,L) | !------- ENDDO | ! OPENMP ENDDO | !------- C----------------------------------------------------------------------- | ! 500 CONTINUE | !$omp parallel do private (IX , JX) C----------------------------------------------------------------------- | ! WP=0.075 | DO 450 K=1,LM C | ! KN=0 | DO J=1,JHL KP=0 | IF (J >= MY_JS_GLB-JBPAD2 .AND. J <= MY_JE_GLB+JTPAD2) THEN DWDTMX=0. | JX = J - MY_JS_GLB + 1 DWDTMN=0. | DO I=1,IM C----------------------------------------------------------------------- | IF (I >= MY_IS_GLB-ILPAD2 .AND. I <= MY_IE_GLB+IRPAD2) THEN !$omp parallel do | IX = I - MY_IS_GLB + 1 !$omp& private (dwdtmn,dwdtmx,dwdtt,hm,kn,kp) | DWDT(IX,JX,K) = MAX(DWDT(IX,JX,k) * HBM3(IX,JX), -0.03) DO 525 L=1,LM | DWDT(IX,JX,K) = MIN(DWDT(IX,JX,K) * HBM3(IX,JX), 0.03) C----------------------------------------------------------------------- | END IF DO J=MYJS,MYJE | END DO DO I=MYIS,MYIE | END IF HM=HTM(I,J,L)*HBM2(I,J) | END DO DWDTT=DWDT(I,J,L)*HM | ! C | DO J=JHH,JM DWDTMX=AMAX1(DWDTT,DWDTMX) | IF (J >= MY_JS_GLB-JBPAD2 .AND. J <= MY_JE_GLB+JTPAD2) THEN DWDTMN=AMIN1(DWDTT,DWDTMN) | JX = J - MY_JS_GLB + 1 C | DO I=1,IM IF(DWDTT.LT.EPSN)THEN | IF (I >= MY_IS_GLB-ILPAD2 .AND. I <= MY_IE_GLB+IRPAD2) THEN DWDTT=EPSN | IX = I - MY_IS_GLB + 1 KN=KN+1 | DWDT(IX,JX,K) = MAX(DWDT(IX,JX,K) * HBM3(IX,JX), -0.03) ENDIF | DWDT(IX,JX,K) = MIN(DWDT(IX,JX,K) * HBM3(IX,JX), 0.03) C | END IF IF(DWDTT.GT.EPSP)THEN | END DO DWDTT=EPSP | END IF KP=KP+1 | END DO ENDIF | ! C | DO J=1,JM DWDT(I,J,L)=(DWDTT/G+1.)*(1.-WP)+PDWDT(I,J,L)*WP | IF (J >= MY_JS_GLB-JBPAD2 .AND. J <= MY_JE_GLB+JTPAD2) THEN ENDDO | JX = J - MY_JS_GLB + 1 ENDDO | DO I=1,IHL C----------------------------------------------------------------------- | IF (I >= MY_IS_GLB-ILPAD2 .AND. I <= MY_IE_GLB+IRPAD2) THEN 525 CONTINUE | IX = I - MY_IS_GLB + 1 C----------------------------------------------------------------------- | DWDT(IX,JX,L) = MAX(DWDT(IX,JX,L) * HBM3(IX,JX), -0.03) GDT=G*DT | DWDT(IX,JX,L) = MIN(DWDT(IX,JX,L) * HBM3(IX,JX), 0.03) GDT2=GDT*GDT | END IF WGHT=0.35 | END DO FFC=-4.*R/GDT2 | END IF C | END DO C----------------------------------------------------------------------- | ! !$omp parallel do | DO J=1,JM !$omp& private (chi,chin,chmod,clim,coff,delp,dfrc,dp,dpstr, | IF (J >= MY_JS_GLB-JBPAD2 .AND. J <= MY_JE_GLB+JTPAD2) THEN !$omp& dptl,dptu,dwdtt,hbm2ij,imax,inot,iter,itmx,jmax, | JX = J - MY_JS_GLB + 1 !$omp& l,lmp,pdp,pnp1,pone,pp1,pstr,rdp,resdl,wil,wres) | DO I=IHH,IM C----------------------------------------------------------------------- | IF (I >= MY_IS_GLB-ILPAD2 .AND. I <= MY_IE_GLB+IRPAD2) THEN DO 600 J=MYJS2,MYJE2 | IX = I - MY_IS_GLB + 1 DO 600 I=MYIS1,MYIE1 | DWDT(IX,JX,K) = MAX(DWDT(IX,JX,K) * HBM3(IX,JX), -0.03) C----------------------------------------------------------------------- | DWDT(IX,JX,K) = MIN(DWDT(IX,JX,K) * HBM3(IX,JX), 0.03) C----------------------------------------------------------------------- | END IF C | END DO LMP=LMH(I,J) | END IF PDP=PDSL(I,J) | END DO RPD=1./PDP | ! C----------------------------------------------------------------------- | 450 END DO PONE(1)=PT | ! PSTR(1)=PT | END IF PNP1(1)=PT | ! CHI(1)=0. | CALL EXCH(DWDT,LM,1,1) C | ! DO L=2,LMP+1 | WA = 0.15 CHI(L)=0. | !------- DPPL=DETA(L-1)*PDP | ! OPENMP RDPP(L-1)=1./DPPL | !------- PONE(L)=PINT(I,J,L) | ! DPSTR=DWDT(I,J,L-1)*DPPL | !$omp parallel do private (ANE , ASE) PSTR(L)=PSTR(L-1)+DPSTR | ! PP1=PNP1(L-1)+DPSTR | DO 500 K=1,LM PNP1(L)=(PP1-PONE(L))*WGHT+PONE(L) | ! TFC=Q(I,J,L-1)*0.608+1. | DO J=MYJS1_P1,MYJE2_P1 TTFC=-CAPA*TFC+1. | DO I=MYIS_P1,MYIE1_P1 COFF(L-1)=T(I,J,L-1)*TTFC*TFC*DPPL*FFC | ANE(I,J) = (DWDT(I+IHE(J),J+1,K) - DWDT(I,J,K)) * HTM(I,J,K) & /((PNP1(L-1)+PNP1(L))*(PNP1(L-1)+PNP1(L))) | & * HTM(I+IHE(J),J+1,K) ENDDO | END DO C | END DO PSTRUP=-(PSTR(1)+PSTR(2)-PONE(1)-PONE(2))*COFF(1) | ! C | DO J=MYJS2_P1,MYJE1_P1 DO L=2,LMP | DO I=MYIS_P1,MYIE1_P1 RDPDN=RDPP(L) | ASE(I,J) = (DWDT(I+IHE(J),J-1,K) - DWDT(I,J,K)) * HTM(I+IHE(J),J-1,K) RDPUP=RDPP(L-1) | & * HTM(I,J,K) C | END DO PSTRDN=-(PSTR(L)+PSTR(L+1)-PONE(L)-PONE(L+1))*COFF(L) | END DO C | ! B1(L)=COFF(L-1)+RDPUP | DO J=MYJS2,MYJE2 B2(L)=(COFF(L-1)+COFF(L))-(RDPUP+RDPDN) | DO I=MYIS1,MYIE1 B3(L)=COFF(L)+RDPDN | DWDT(I,J,K) = (ANE(I,J) - ANE(I+IHW(J),J-1) + ASE(I,J) - ASE(I+IHW(J),J+1)) C0(L)=PSTRUP+PSTRDN | & * WA * HBM2(I,J) + DWDT(I,J,K) C | END DO PSTRUP=PSTRDN | END DO ENDDO | ! C | 500 END DO B1(2)=0. | ! B2(LMP)=B2(LMP)+B3(LMP) | WP = 0.075 C---------------ELIMINATION--------------------------------------------- | KN = 0 DO L=3,LMP | KP = 0 TMP=-B1(L)/B2(L-1) | DWDTMX = 0. B2(L)=B3(L-1)*TMP+B2(L) | DWDTMN = 0. C0(L)=C0(L-1)*TMP+C0(L) | !------- ENDDO | ! OPENMP C | !------- CHI(1)=0. | ! C---------------BACK SUBSTITUTION-------------------------------------- | !$omp parallel do CHI(LMP)=C0(LMP)/B2(LMP) | !$omp private (DWDTMN , DWDTMX , DWDTT , HM , KN , KP) CHI(LMP+1)=CHI(LMP) | ! C | DO 525 K=1,LM DO L=LMP-1,2,-1 | ! CHI(L)=(-B3(L)*CHI(L+1)+C0(L))/B2(L) | DO J=MYJS,MYJE ENDDO | DO I=MYIS,MYIE C----------------------------------------------------------------------- | HM = HTM(I,J,K) * HBM2(I,J) HBM2IJ=HBM2(I,J) | DWDTT = DWDT(I,J,K) * HM DPTU=0. | ! FCT=0.5/CP*HBM2IJ | DWDTMX = AMAX1(DWDTT,DWDTMX) C | DWDTMN = AMIN1(DWDTT,DWDTMN) DO L=1,LMP | ! DPTL=(CHI(L+1)+PSTR(L+1)-PINT(I,J,L+1))*HBM2IJ | IF (DWDTT < EPSN) THEN PINT(I,J,L+1)=PINT(I,J,L+1)+DPTL | DWDTT = EPSN T(I,J,L)=(DPTU+DPTL)*RTOP(I,J,L)*FCT+T(I,J,L) | KN = KN + 1 DELP=(PINT(I,J,L+1)-PINT(I,J,L))*RDPP(L) | END IF W(I,J,L)=((DELP-DWDT(I,J,L))*GDT+W(I,J,L))*HBM2IJ | ! DWDT(I,J,L)=(DELP-1.)*HBM2IJ+1. | IF (DWDTT > EPSP) THEN C | DWDTT = EPSP DPTU=DPTL | KP = KP + 1 ENDDO | END IF C | ! DO L=LMP+1,LM | DWDT(I,J,K) = (DWDTT/G+1.) * (1.-WP) + PDWDT(I,J,K) * WP PINT(I,J,L+1)=PDSL(I,J)*DETA(L)+PINT(I,J,L) | END DO ENDDO | END DO C----------------------------------------------------------------------- | ! 600 CONTINUE | 525 END DO C----------------------------------------------------------------------- | ! C write(6,*) 'for PE: ', MYPE, 'DWDT extrema: ', DWDTMN,DWDTMX | GDT = G * DT RETURN | GDT2 = GDT * GDT END | WGHT = 0.35 > FFC = -4. * R / GDT2 > !------- > ! OPENMP > !------- > ! > !$omp parallel do > !$omp private (CHI , CHIN , CHMOD , CLIM , COFF , DELP , DFRC , DP , > !$omp DPSTR , DPTL , DPTU , DWDTT , HBM2IJ , IMAX , INOT , ITER , > !$omp ITMX , JMAX , K , LMP , PDP , PNP1 , PONE , PP1 , > !$omp PSTR , RDP , RESDL , WIL , WRES ) > ! > DO 600 J=MYJS2,MYJE2 > DO 600 I=MYIS1,MYIE1 > ! > LMP = LMH(I,J) > PDP = PDSL(I,J) > RPD = 1. / PDP > ! > PONE(1) = PT > PSTR(1) = PT > PNP1(1) = PT > CHI(1) = 0. > ! > DO K=2,LMP+1 > CHI(K) = 0. > DPPL = DETA(K-1) * PDP > RDPP(K-1) = 1. / DPPL > PONE(K) = PINT(I,J,K) > DPSTR = DWDT(I,J,K-1) * DPPL > PSTR(K) = PSTR(K-1) + DPSTR > PP1 = PNP1(K-1) + DPSTR > PNP1(K) = (PP1-PONE(K)) * WGHT + PONE(K) > TFC = Q(I,J,K-1) * 0.608 + 1. > TTFC = -CAPA * TFC + 1. > COFF(K-1) = T(I,J,K-1) * TTFC * TFC * DPPL * FFC / ((PNP1(K-1)+PNP1(K)) * > & (PNP1(K-1)+PNP1(K))) > END DO > ! > PSTRUP = -(PSTR(1) + PSTR(2) - PONE(1) - PONE(2)) * COFF(1) > ! > DO K=2,LMP > RDPDN = RDPP(K) > RDPUP = RDPP(K-1) > ! > PSTRDN = -(PSTR(K) + PSTR(K+1) - PONE(K) - PONE(K+1)) * COFF(K) > ! > B1(K) = COFF(K-1) + RDPUP > B2(K) = (COFF(K-1) + COFF(K)) - (RDPUP+RDPDN) > B3(K) = COFF(K) + RDPDN > C0(K) = PSTRUP + PSTRDN > ! > PSTRUP = PSTRDN > END DO > ! > B1(2) = 0. > B2(LMP) = B2(LMP) + B3(LMP) > !------------ > ! ELIMINATION > !------------ > DO K=3,LMP > TMP = -B1(K) / B2(K-1) > B2(K) = B3(K-1) * TMP + B2(K) > C0(K) = C0(K-1) * TMP + C0(K) > END DO > ! > CHI(1) = 0. > !------------------ > ! BACK SUBSTITUTION > !------------------ > CHI(LMP) = C0(LMP) / B2(LMP) > CHI(LMP+1) = CHI(LMP) > ! > DO K=LMP-1,2,-1 > CHI(K) = (-B3(k) * CHI(K+1) + C0(K)) / B2(K) > END DO > ! > HBM2IJ = HBM2(I,J) > DPTU = 0. > FCT = 0.5 / CP * HBM2IJ > ! > DO K=1,LMP > DPTL = (CHI(K+1) + PSTR(K+1) - PINT(I,J,K+1)) * HBM2IJ > PINT(I,J,K+1) = PINT(I,J,K+1) + DPTL > T (I,J,K) = (DPTU+DPTL) * RTOP(I,J,K) * FCT + T(I,J,K) > DELP = (PINT(I,J,K+1) - PINT(I,J,K)) * RDPP(K) > W(I,J,K) = ((DELP - DWDT(I,J,K)) * GDT + W(I,J,K)) * HBM > DWDT(I,J,K) = (DELP-1.) * HBM2IJ + 1. > ! > DPTU = DPTL > END DO > ! > DO K=LMP,LM > PINT(I,J,K+1) = PDSL(I,J) * DETA(K) + PINT(I,J,K) > END DO > ! > 600 END DO > ! > RETURN > ! > END SUBROUTINE EPS