SUBROUTINE COMP implicit none C **** ADVANCE PSI IN TIME include "params.h" include "strt.h" include "blnk.h" include "ccomp.h" include "cons.h" include "index.h" include "buff.h" include "phys.h" include "cmpbnd.h" include "crates.h" ! ! Local: real :: AK0(2,2),PHI(2,3),DELTA(2,2),tau,t00,small integer :: i,l,ll,nps1k,km,kp,npsl,m,ntk,k,kb,nwk,nrjk,npsk, | npjp1k,npsnmk,npjm1k,npjm2k,npsdhk,kk,npsm,npsnpk,npsmnk,len,k1, | k2,n,nps1mk,nps2mk,nps2k,npjp2k ! DATA PHI/0.,0.673,1.35,0.,1.11,0.769/,TAU/1.86E+3/, ADELTA/1.,0.,0.,1./, BT00/273./ DATA SMALL/1.E-20/ DO 34 I=1,LEN1 T1(I)=1. 34 CONTINUE IF(ICOMP.EQ.1)CALL DFACT(T1,J,IMAXP4) C **** CALCULATE EMBAR NPS1K=NJ+NPS NPS2K=NPS1K+KMAXP1 DO 1 I=1,LEN3 EMBAR(I,1)=1./(F(I,NPS1K)/RMASS(1)+F(I,NPS2K)/RMASS(2)+(1.- AF(I,NPS1K)-F(I,NPS2K))/RMASS(3)) 1 CONTINUE C CALCULATE PS0 AND EMBAR0 DO 30 I=1,LEN1 PS0(I,1) = B(I,1,1)*F(I,NPS1K)+B(I,1,2)*F(I,NPS2K)+FB(I,1) PS0(I,2) = B(I,2,1)*F(I,NPS1K)+B(I,2,2)*F(I,NPS2K)+FB(I,2) EMBAR0(I)=1./(PS0(I,1)/RMASS(1)+PS0(I,2)/RMASS(2)+(1.-PS0(I,1) A-PS0(I,2))/RMASS(3)) C **** WKS4=.5*(DMBAR/DZ)/MBAR WKS4(I)=(EMBAR(I,1)-EMBAR0(I))/(C(3)*(EMBAR0(I)+EMBAR(I,1))) 30 CONTINUE C **** CALCULTE EP AND AK AT LEVEL 1/2, SET Z(0) TO ZERO KM=1 KP=2 C **** EP(1/2) DO 2 L=1,2 DO 2 I=1,LEN1 EP(I,L,KP)=1.-(2./(EMBAR0(I)+EMBAR(I,1)))*(RMASS(L)+(EMBAR(I,1)- AEMBAR0(I))/C(3)) Z(I,1,L)=0. 2 CONTINUE DO 3 L=1,2 LL=3-L NPSL=NJ+NPS+(L-1)*KMAXP1 DO 3 M=1,2 DO 3 I=1,LEN1 AK(I,L,M,KP)=-DELTA(L,M)*(PHI(LL,3)+(PHI(LL,L)-PHI(LL,3))*.5* A(PS0(I,L)+F(I,NPSL)))-(1.-DELTA(L,M))*(PHI(L,M)-PHI(L,3))* B.5*(PS0(I,L)+F(I,NPSL)) 3 CONTINUE C **** WKS1=MBAR/M3*(T00/(T0+T))*0.25/(TAU*DET) AT LEVEL 1/2 NTK=NJ+NT+KMAX DO 4 I=1,LEN1 WKS1(I)=0.5*(EMBAR0(I)+EMBAR(I,1))/RMASS(3)*(T00/(T0(1)+F(I,NTK))) A**0.25/(TAU*(AK(I,1,1,KP)*AK(I,2,2,KP)-AK(I,1,2,KP)*AK(I,2,1,KP))) 4 CONTINUE C **** COMPLETE CALCULATION OF AK(1/2) DO 5 L=1,2 DO 5 M=1,2 DO 5 I=1,LEN1 AK(I,L,M,KP)=AK(I,L,M,KP)*WKS1(I) C **** SET GAMA(0) TO 0. GAMA(I,1,L,M)=0. 5 CONTINUE C **** MAIN HEIGHT DO LOOP DO 6 K=1,KMAX C **** CYCLE K LEVELS KB=KM KM=KP KP=KB C **** EP(K+1/2) DO 7 L=1,2 DO 7 I=1,LEN1 EP(I,L,KP)=1.-(2./(EMBAR(I,K)+EMBAR(I,K+1)))*(RMASS(L)+ A(EMBAR(I,K+1)-EMBAR(I,K))/C(3)) 7 CONTINUE C **** AK(K+1/2) DO 8 L=1,2 LL=3-L NPSL=NJ+NPS+(L-1)*KMAXP1+K DO 8 M=1,2 DO 8 I=1,LEN1 AK(I,L,M,KP)=-DELTA(L,M)*(PHI(LL,3)+(PHI(LL,L)-PHI(LL,3))*.5* A(F(I,NPSL-1)+F(I,NPSL)))-(1.-DELTA(L,M))*(PHI(L,M)-PHI(L,3))* B.5*(F(I,NPSL-1)+F(I,NPSL)) 8 CONTINUE C **** WKS1=MBAR/M3*(T00/(T0+T))**0.25/(TAU*DET(ALFA)) NTK=NJ+NT+K DO 9 I=1,LEN1 WKS1(I)=0.5*(EMBAR(I,K)+EMBAR(I,K+1))/RMASS(3)*(T00/(T0(K+1)+.5* A(F(I,NTK-1)+F(I,NTK))))**0.25/(TAU*(AK(I,1,1,KP)*AK(I,2,2,KP)- BAK(I,1,2,KP)*AK(I,2,1,KP))) C **** EDDY DIFFUSION TERMS IN WKS3 AND WKS4 WKS3(I)=WKS4(I) WKS4(I)=(EMBAR(I,K+1)-EMBAR(I,K))/(C(3)*(EMBAR(I,K)+EMBAR(I,K+1))) 9 CONTINUE C **** FINISH CALCULATING AK(K+1/2) AND GENERATE PK, QK, RK NWK=NJ+NW+K NRJK=NRJ+K DO 10 L=1,2 DO 10 M=1,2 DO 10 I=1,LEN1 AK(I,L,M,KP)=AK(I,L,M,KP)*WKS1(I) PK(I,L,M)=(AK(I,L,M,KM)*(1./C(3)+EP(I,M,KM)/2.)-EXPS(K) A*(C(87)*DIFK(K)*T1(I)*(1./C(3)-WKS3(I))+0.25*(F(I,NWK-1)+ BF(I,NWK)))*DELTA(L,M))/C(3) RK(I,L,M)=(AK(I,L,M,KP)*(1./C(3)-EP(I,M,KP)/2.)-EXPS(K) A*(C(86)*DIFK(K+1)*T1(I)*(1./C(3)+WKS4(I))-0.25*(F(I,NWK-1)+ BF(I,NWK)))*DELTA(L,M))/C(3) QK(I,L,M)=-(AK(I,L,M,KM)*(1./C(3)-EP(I,M,KM)/2.)+AK(I,L,M,KP)* A(1./C(3)+EP(I,M,KP)/2.))/C(3)+EXPS(K)*(((C(86)*DIFK(K+1)* B(1./C(3)-WKS4(I))+C(87)*DIFK(K)*(1./C(3)+WKS3(I)))*T1(I)/C(3)+ BC(7))*DELTA(L,M)-FS(I,K,L,M)) 10 CONTINUE C **** CALCULATE HORIZONTAL ADVECTION IN FK ARRAY DO 11 L=1,2 NPSL=NPS+(L-1)*KMAXP1+K-1 CALL ADVECL(NPSL,FK(1,L),K) 11 CONTINUE C **** C **** ADD EXPLICIT SOURCE TERMS TO FK C **** DO 211 I = 1,LEN1 FK(I,1) = FK(I,1)-FS(I,K,1,0) FK(I,2) = FK(I,2)-FS(I,K,2,0) 211 CONTINUE C **** SHAPIRO SMOOTHER DO 111 L=1,2 NPSK=NPSNM+(L-1)*KMAXP1+K-1 NPJP2K=NJP2+NPSK NPJP1K=NJP1+NPSK NPSNMK=NJ+NPSK NPJM1K=NJM1+NPSK NPJM2K=NJM2+NPSK DO 111 I=1,LEN1 WKV1(I,L)=F(I,NPSNMK)-C(26)*(F(I,NPJP2K)+F(I,NPJM2K)-4.* A(F(I,NPJP1K)+F(I,NPJM1K))+6.*F(I,NPSNMK)) 111 CONTINUE DO 112 L=1,2 DO 112 I=3,IMAXP2 WKV2(I,L)=WKV1(I,L)-C(26)*(WKV1(I+2,L)+WKV1(I-2,L)-4.* A(WKV1(I+1,L)+WKV1(I-1,L))+6.*WKV1(I,L)) 112 CONTINUE C **** COMPLETE CALCULATION OF RHS IN FK DO 12 L=1,2 ! NPSDHK=NPSDH+(L-1)*KMAX+K-1 NPSDHK=NPSDH+(L-1)*KMAXP1+K-1 DO 12 I=3,IMAXP2 FK(I,L)=EXPS(K)*(WKV2(I,L)*C(7)-FK(I,L)+F(I,NPSDHK)) 12 CONTINUE C **** INSERT PERIODIC POINTS DO 121 L=1,2 DO 121 I=1,2 FK(I,L)=FK(I+IMAX,L) FK(I+IMAXP2,L)=FK(I+2,L) 121 CONTINUE IF(K.EQ.1) GO TO 13 IF(K.EQ.KMAX) GO TO 14 GO TO 15 13 CONTINUE C **** LOWER BNDRY DO 16 L=1,2 DO 16 M=1,2 DO 16 KK=1,2 DO 16 I=1,LEN1 QK(I,L,M) = QK(I,L,M)+PK(I,L,KK)*B(I,KK,M) 16 CONTINUE DO 33 L=1,2 DO 33 M=1,2 DO 33 I=1,LEN1 FK(I,L) = FK(I,L)-PK(I,L,M)*FB(I,M) PK(I,L,M)=0. 33 CONTINUE GO TO 15 14 CONTINUE C **** UPPER BNDRY DO 17 L=1,2 DO 17 M=1,2 DO 17 I=1,LEN1 QK(I,L,M)=QK(I,L,M)+(1.+.5*EP(I,M,KP)*C(3))/(1.-.5*EP(I,M,KP)* AC(3))*RK(I,L,M) RK(I,L,M)=0. 17 CONTINUE 15 CONTINUE C **** QK=ALFAK=QK-PK*GAMA(K-1) DO 18 L=1,2 DO 18 M=1,2 DO 18 KK=1,2 DO 18 I=1,LEN1 QK(I,L,M)=QK(I,L,M)-PK(I,L,KK)*GAMA(I,K,KK,M) 18 CONTINUE C **** WKS1=DET(ALFA) DO 19 I=1,LEN1 WKS1(I)=QK(I,1,1)*QK(I,2,2)-QK(I,1,2)*QK(I,2,1) 19 CONTINUE C **** WKM1=ALFAI DO 20 L=1,2 LL=3-L DO 20 M=1,2 DO 20 I=1,LEN1 WKM1(I,L,M)=(DELTA(L,M)*QK(I,LL,LL)-(1.-DELTA(L,M))*QK(I,L,M))/ AWKS1(I) 20 CONTINUE C **** GAMA(K+1)=ALFAI*RK C **** WKV1=FK-PK*Z(K) DO 21 L=1,2 DO 22 I=1,LEN1 WKV1(I,L)=FK(I,L) 22 CONTINUE DO 21 M=1,2 DO 23 I=1,LEN1 GAMA(I,K+1,L,M)=0. WKV1(I,L)=WKV1(I,L)-PK(I,L,M)*Z(I,K,M) 23 CONTINUE DO 21 KK=1,2 DO 21 I=1,LEN1 GAMA(I,K+1,L,M)=GAMA(I,K+1,L,M)+WKM1(I,L,KK)*RK(I,KK,M) 21 CONTINUE C **** Z(K+1)=WKM1*WKV1 DO 231 L=1,2 DO 232 I=1,LEN1 Z(I,K+1,L)=0. 232 CONTINUE DO 231 M=1,2 DO 231 I=1,LEN1 Z(I,K+1,L)=Z(I,K+1,L)+WKM1(I,L,M)*WKV1(I,M) 231 CONTINUE 6 CONTINUE C **** SET PSNP(KMAXP1) TO ZERO DO 24 L=1,2 NPSL=NJNP+NPS+(L-1)*KMAXP1+KMAX DO 24 I=1,LEN1 F(I,NPSL)=0. 24 CONTINUE C **** DOWNWARD SWEEP DO 25 KK=1,KMAX K=KMAX+1-KK DO 25 L=1,2 NPSL=NJNP+NPS+(L-1)*KMAXP1+(K-1) DO 26 I=1,LEN1 F(I,NPSL)=Z(I,K+1,L) 26 CONTINUE DO 25 M=1,2 NPSM=NJNP+NPS+(M-1)*KMAXP1+K DO 25 I=1,LEN1 F(I,NPSL)=F(I,NPSL)-GAMA(I,K+1,L,M)*F(I,NPSM) 25 CONTINUE C **** INSERT VALUE OF PS(KMAXP1) USING BNDRY CONDITION DO 27 L=1,2 NPSL=NJNP+NPS+(L-1)*KMAXP1+KMAX DO 27 I=1,LEN1 F(I,NPSL)=(1.+.5*EP(I,L,KP)*C(3))/(1.-.5*EP(I,L,KP)*C(3))* AF(I,NPSL-1) 27 CONTINUE C **** FILTER PSI IF(KUT(J).GE.IMAXH) GO TO 28 NPSK=NJNP+NPS CALL FILTER(NPSK,2*KMAXP1,KUT(J)) 28 CONTINUE C **** TIME SMOOTHING OF PSI NPSNMK=NJ+NPSNM NPSK=NJ+NPS NPSNPK=NJNP+NPS NPSMNK=NJNP+NPSNM LEN=2.*LEN3 DO 29 I=1,LEN F(I,NPSMNK)=C(30)*F(I,NPSK)+C(31)*(F(I,NPSNMK)+F(I,NPSNPK)) 29 CONTINUE C **** INSERT PERIODIC POINTS K1=NJNP+NPS K2=K1+2*KMAXP1-1 DO 31 N=1,2 DO 32 I=1,2 DO 32 K=K1,K2 F(I,K)=F(I+IMAX,K) F(I+IMAXP2,K)=F(I+2,K) 32 CONTINUE K1=NJNP+NPSNM K2=K1+2*KMAXP1-1 31 CONTINUE C **** INSURE NON-NEGATIVE PSI NPS1K=NJNP+NPS NPS2K=NJNP+NPS2 NPS1MK=NJNP+NPSNM NPS2MK=NJNP+NPS2NM DO 41 I=1,LEN3 F(I,NPS1K)=merge(F(I,NPS1K),SMALL,F(I,NPS1K)-SMALL>=0.) F(I,NPS2K)=merge(F(I,NPS2K),SMALL,F(I,NPS2K)-SMALL>=0.) F(I,NPS1MK)=merge(F(I,NPS1MK),SMALL,F(I,NPS1MK)-SMALL>=0.) F(I,NPS2MK)=merge(F(I,NPS2MK),SMALL,F(I,NPS2MK)-SMALL>=0.) F(I,NPS1K)=F(I,NPS1K)*merge(1.,(1.-SMALL)/(F(I,NPS1K)+F(I,NPS2K)), 11.-SMALL-F(I,NPS1K)-F(I,NPS2K)>=0.) F(I,NPS2K)=F(I,NPS2K)*merge(1.,(1.-SMALL)/(F(I,NPS1K)+F(I,NPS2K)), 11.-SMALL-F(I,NPS1K)-F(I,NPS2K)>=0.) F(I,NPS1MK)=F(I,NPS1MK)*merge(1.,(1.-SMALL)/(F(I,NPS1MK)+ 1F(I,NPS2MK)),1.-SMALL-F(I,NPS1MK)-F(I,NPS2MK)>=0.) F(I,NPS2MK)=F(I,NPS2MK)*merge(1.,(1.-SMALL)/(F(I,NPS1MK)+ 1F(I,NPS2MK)),1.-SMALL-F(I,NPS1MK)-F(I,NPS2MK)>=0.) 41 CONTINUE RETURN END C