subroutine slp implicit none include 'param_o.h' include 'vrbls.comm' include 'extra.comm' include 'mapot.comm' include 'const.h' include 'dynam.comm' include 'masks.comm' include 'params' integer,parameter::nrlx1=500 integer::i,j,n,l,lhmnt,lmap1,ll,nrlx,kmn,kmntm,lmst,li,kmm, & km,n1 integer::isCorner2,lma real,dimension(0:im+1,0:jm+1,nm)::ttv,pdsl1,pbi,slpx integer,dimension((im+2)*(jm+2)*nm)::imnt,jmnt,nmnt real::tgss,pbin,phbi,ptin,trtv,phti,dposp,prbin,prtin,alpp1,slop, & slpp,ttt real,parameter::OVERRC=1.50,AD05=OVERRC*0.05,CFT0=OVERRC-1., & r=287.04,rog=r/9.8 logical:: stdrd if(sigma)then stdrd=.true. else stdrd=.false. endif ! !----------------------------------------------------------------------- ! !*** FIND THE MOST ELEVATED GLOBAL LAYER WHERE THERE IS LAND ! !----------------------------------------------------------------------- DO L=1,LM do n=1,nm DO J=1,jm DO I=1,IM IF(htm(I,J,n,l).LT.0.5)GO TO 666 ENDDO ENDDO ENDDO ! !*** IF WE GET TO HERE, WE HAVE ALL ATMOSPHERE ! GOTO 667 666 CONTINUE ! !*** IF WE GET HERE, WE HAVE FOUND A NON_ATM POINT ! LHMNT=L GOTO 668 667 CONTINUE ENDDO ! LHMNT=LM+1 !cccccc GO TO 669 668 CONTINUE !----------------------------------------------------------------------- !*** !*** INITIALIZE ARRAYS. LOAD SLP ARRAY WITH SURFACE PRESSURE. !*** do n=1,nm DO J=1,jm DO I=1,IM PSLP(I,J,n)=0. TTV(I,J,n)=0. ENDDO ENDDO ENDDO ! do 110 n=1,nm DO 110 J=1,jm DO 110 I=1,IM PDSL1(I,J,n)=RES(I,J,n)*PD(I,J,n) PSLP(I,J,n)=PD(I,J,n)+PT PBI (I,J,n)=PSLP(I,J,n) 110 CONTINUE if(.not.stdrd)then LL=LM ! do n=1,nm DO J=1,jm DO I=1,im IF(HTM(I,J,n,LL).LE.0.5) THEN TGSS=FIS(I,J,n)/(R*ALOG((PDSL1(I,J,n)+PT)/(PD(I,J,n)+PT))) LMAP1=LMH(I,J,n)+1 ! if(n.eq.3.and.i.eq.30.and.j.eq.101)then ! print *,i,j,tgss,fis(i,j,n),r,pdsl1(i,j,n),pd(i,j,n) ! stop endif DO 260 L=LMAP1,LM T(I,J,n,L)=TGSS 260 CONTINUE ! ENDIF ENDDO ENDDO ENDDO !---------------------------------------------------------------- ! !*** CREATE A TEMPORARY TV ARRAY, AND FOLLOW BY SEQUENTIAL !*** OVERRELAXATION, DOING NRLX PASSES. ! !---------------------------------------------------------------- NRLX=NRLX1 ! !---------------------------------------------------------------- !---------------------------------------------------------------- DO 300 L=LHMNT,LM !---------------------------------------------------------------- !---------------------------------------------------------------- ! KMN=0 KMNTM=0 ! do 240 n=1,nm DO 240 J=1,jm DO 240 I=1,im IF(HTM(I,J,n,L).GT.0.5)GO TO 240 KMN=KMN+1 IMNT(KMN)=I JMNT(KMN)=J nmnt(kmn)=n 240 CONTINUE ! KMNTM=KMN ! do 270 n=1,nm DO 270 J=1,jm DO 270 I=1,IM TTV(I,J,n)=T(I,J,n,L) 270 CONTINUE ! !---------------------------------------------------------------- !*** FOR GRID BOXES NEXT TO MOUNTAINS REPLACE TTV BY AN "EQUIVALENT" !*** TV, ONE WHICH CORRESPONDS TO THE CHANGE IN P BETWEEN REFERENCE !*** INTERFACE GEOPOTENTIALS, INSTEAD OF BETWEEN LAYER INTERFACES !---------------------------------------------------------------- ! do n=1,nm DO J=1,jm DO I=1,im IF(HTM(I,J,n,L).GT.0.5.AND. & HTM(I,J-1,n,L)*HTM(I+1,J-1,n,L) & *HTM(I,J+1,n,L)*HTM(I+1,J+1,n,L) & *HTM(I-1 ,J ,n,L)*HTM(I+1 ,J ,n,L) & *HTM(I-1 ,J-1,n,L)*HTM(I-1 ,J+1,n,L).LT.0.5)THEN LMST=LMH(I,J,n) !*** !*** FIND P AT THE REFERENCE INTERFACE GEOPOTENTIAL AT THE BOTTOM !*** PBIN=PT+PD(I,J,n) PHBI=DFL(LMST+1) ! DO LI=LMST,1,-1 PTIN=PBIN-DETA(LI)*PD(I,J,n)*RES(I,J,n) TRTV=2.*R*T(I,J,n,LI)*(1.+0.608*Q(I,J,n,LI)) PHTI=PHBI+TRTV*(PBIN-PTIN)/(PBIN+PTIN) if(i.eq.17.and.j.eq.89.and.n.eq.5)then ! print *,l,pbin,ptin,deta(li),res(i,j,n),pd(i,j,n) endif IF(PHTI.GE.DFL(L+1))GO TO 273 PBIN=PTIN PHBI=PHTI ENDDO ! 273 DPOSP=(PHTI-DFL(L+1))/TRTV PRBIN=(1.+DPOSP)/(1.-DPOSP)*PTIN !*** !*** FIND P AT THE REFERENCE INTERFACE GEOPOTENTIAL AT THE TOP !*** PBIN=PT+PD(I,J,n) PHBI=DFL(LMST+1) ! DO LI=LMST,1,-1 PTIN=PBIN-DETA(LI)*PD(I,J,n)*RES(I,J,n) TRTV=2.*R*T(I,J,n,LI)*(1.+0.608*Q(I,J,n,LI)) PHTI=PHBI+TRTV*(PBIN-PTIN)/(PBIN+PTIN) IF(PHTI.GE.DFL(L))GO TO 275 PBIN=PTIN PHBI=PHTI ENDDO ! 275 DPOSP=(PHTI-DFL(L))/TRTV PRTIN=(1.+DPOSP)/(1.-DPOSP)*PTIN ! TTV(I,J,n)=(DFL(L)-DFL(L+1))/(2.*R)*(PRBIN+PRTIN)/(PRBIN-PRTIN) ENDIF ENDDO ENDDO ENDDO !---------------------------------------------------------------- KMM=KMNTM !---------------------------------------------------------------- !*** !*** HERE IS THE RELAXATION LOOP !*** !---------------------------------------------------------------- DO 285 N1=1,NRLX ! CALL bocoh(TTV,im,jm,nm) ! Exchange haloes ! DO 280 KM=1,KMM I=IMNT(KM) J=JMNT(KM) n=nMNT(KM) if(isCorner2(i,j,n,im,jm,nm).eq.0)then TTV(I,J,n)=AD05*(4.*(TTV(I,J-1,n)+TTV(I+1,J,n) & +TTV(I,J+1,n)+TTV(I-1,J,n)) & +TTV(I-1,J-1,n) +TTV(I-1,J+1,n) & +TTV(I+1,J-1,n) +TTV(I+1,J+1,n)) & -CFT0*TTV(I,J,n) else TTV(i,j,n)=0.25*(TTV(I,J-1,n)+TTV(I+1,J,n) & +TTV(I,J+1,n)+TTV(I-1,J,n)) endif 280 CONTINUE ! 285 CONTINUE !---------------------------------------------------------------- ! if(l.eq.36)then ! print *,'t(30,101,3),',t(30,101,3,36),ttv(30,101,3) endif DO 290 KM=1,KMM I=IMNT(KM) J=JMNT(KM) n=nMNT(KM) if(l.eq.36)then ! print *,i,j,n,km endif T(I,J,n,L)=TTV(I,J,n) 290 CONTINUE ! 300 CONTINUE ! VALUES FOR IMNT AND JMNT ARE FOR LAYER LM - THIS IS WHAT WE WANT KMM=KMNTM ! DO 320 KM=1,KMM I=IMNT(KM) J=JMNT(KM) n=nMNT(KM) LMAP1=LMH(I,J,n)+1 PBIN=PT+PD(I,J,n) ! DO L=LMAP1,LM PTIN=PBIN DPOSP=(DFL(L)-DFL(L+1))/(2.*R*T(I,J,n,L)) PBIN=(1.+DPOSP)/(1.-DPOSP)*PTIN if((i.eq.35.or.i.eq.36).and.j.eq.58.and.n.eq.3)then ! print *,i,j,n,l,pd(i,j,n),ptin,dposp ! stop endif ENDDO ! PSLP(I,J,n)=PBIN if(i.eq.70.and.j.eq.186.and.n.eq.4)then ! print *,i,j,n,lmap1,lm ! stop endif if(pslp(i,j,n).gt.123000)then ! print *,i,j,n,pslp(i,j,n) ! stop endif 320 CONTINUE ! stop else write(6,*) 'doing standard reduction!!!' DO 410 n=1,nM DO 410 J=1,JM DO 410 I=1,IM IF(FIS(I,J,n).GE.1.)THEN LMA=LMH(I,J,n) ALPP1=ALOG(PDSL1(I,J,n)*ETA(LMA+1)+PT) SLOP=0.0065*ROG*T(I,J,n,LMA) IF(SLOP.LT.0.50)THEN SLPP=ALPP1+FIS(I,J,n)/(R*T(I,J,n,LMA)) ELSE TTT=-(ALOG(PDSL1(I,J,n)*ETA(LMA)+PT)+ALPP1) & *SLOP*0.50+T(I,J,n,LMA) SLPP=(-TTT+SQRT(TTT*TTT+2.*SLOP* & (FIS(I,J,n)/R+ & (TTT+0.50*SLOP*ALPP1)*ALPP1)))/SLOP ENDIF PSLP(I,J,n)=EXP(SLPP) ENDIF 410 CONTINUE endif do 440 n=1,nm DO 440 J=1,JM DO 440 I=1,IM SLPX(I,J,n)=PSLP(I,J,n) 440 CONTINUE !C !C !C!$OMP parallel do private(ihh2) do 460 n=1,nm DO 460 J=2,JM-1 DO 460 I=2,im-1 !C !C*** EXTRA AVERAGING UNDER MOUNTAINS TAKEN OUT, FM, MARCH 96 !C SLPX(I,J,n)=0.125*(PSLP(I,J-1,n)+PSLP(I,J+1,n) & +PSLP(I+1,J,n)+PSLP(I-1,J,n) & +4.*PSLP(I,J,n)) 460 CONTINUE !C !C!$OMP parallel do do n=1,nm DO J=1,JM DO I=1,IM PSLP(I,J,n)=SLPX(I,J,n) ENDDO ENDDO ENDDO end subroutine slp