#include "dims.h" SUBROUTINE COMP use input_module,only: difhor use cons_module,only: len1,len3,imax,imaxh,imaxp2,imaxp4,kmax, | kmaxp1,difk,expz,t0,kut,rmass,rmassinv,shapiro,dtsmooth, | dtsmooth_div2,expzmid_inv,expzmid,dtx2inv implicit none C **** ADVANCE PSI IN TIME #include "params.h" #include "fcom.h" #include "ccomp.h" #include "index.h" #include "buff.h" #include "phys.h" #include "cmpbnd.h" #include "crates_tdep.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 PHI(:,1)=(/0.,0.673/) PHI(:,2)=(/1.35,0./) PHI(:,3)=(/1.11,0.769/) TAU=1.86e+3 DELTA(:,1)=(/1.,0./) DELTA(:,2)=(/0.,1./) T00=273. ! SMALL = 1.e-20 ! SMALL = 1.e-6 SMALL = 1.e-10 ! DO 34 I=1,LEN1 T1(I)=1. 34 CONTINUE IF (difhor.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)*rmassinv(1)+F(I,NPS2K)*rmassinv(2)+ | (1.-F(I,NPS1K)-F(I,NPS2K))*rmassinv(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)*rmassinv(1)+PS0(I,2)*rmassinv(2)+ | (1.-PS0(I,1)-PS0(I,2))*rmassinv(3)) C **** WKS4=.5*(DMBAR/DZ)/MBAR WKS4(I)=(EMBAR(I,1)-EMBAR0(I))/(dz*(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))/dz) 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))*rmassinv(3)*(T00/(T0(1)+ | F(I,NTK)))**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))/dz) 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))*rmassinv(3)*(T00/(T0(K+1)+ | .5*(F(I,NTK-1)+F(I,NTK))))**0.25/(TAU*(AK(I,1,1,KP)* | AK(I,2,2,KP)-AK(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))/(dz*(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./dz+EP(I,M,KM)/2.)-expz(K) A*(expzmid_inv*DIFK(K)*T1(I)*(1./dz-WKS3(I))+0.25*(F(I,NWK-1)+ BF(I,NWK)))*DELTA(L,M))/dz RK(I,L,M)=(AK(I,L,M,KP)*(1./dz-EP(I,M,KP)/2.)-expz(K) A*(expzmid*DIFK(K+1)*T1(I)*(1./dz+WKS4(I))-0.25*(F(I,NWK-1)+ BF(I,NWK)))*DELTA(L,M))/dz QK(I,L,M)=-(AK(I,L,M,KM)*(1./dz-EP(I,M,KM)/2.)+AK(I,L,M,KP)* A(1./dz+EP(I,M,KP)/2.))/dz+expz(K)*(((expzmid*DIFK(K+1)* B(1./dz-WKS4(I))+expzmid_inv*DIFK(K)*(1./dz+WKS3(I)))*T1(I)/dz+ Bdtx2inv)*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)-shapiro*(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)-shapiro*(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)=expz(K)*(WKV2(I,L)*dtx2inv-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)*dz)/(1.-.5*EP(I,M,KP)* Adz)*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)*dz)/(1.-.5*EP(I,L,KP)*dz)* 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)=dtsmooth*F(I,NPSK)+dtsmooth_div2*(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 if (f(i,nps1k) < small) f(i,nps1k) = small if (f(i,nps2k) < small) f(i,nps2k) = small if (f(i,nps1mk) < small) f(i,nps1mk) = small if (f(i,nps2mk) < small) f(i,nps2mk) = small if (1.-SMALL-F(I,NPS1K)-F(I,NPS2K) < 0.) then f(i,nps1k) = f(i,nps1k)*((1.-SMALL)/(F(I,NPS1K)+F(I,NPS2K))) f(i,nps2k) = f(i,nps2k)*((1.-SMALL)/(F(I,NPS1K)+F(I,NPS2K))) endif if (1.-SMALL-F(I,NPS1MK)-F(I,NPS2MK) < 0.) then f(i,nps1mk) = f(i,nps1mk)*((1.-SMALL)/(F(I,NPS1MK)+ | F(I,NPS2MK))) f(i,nps2mk) = f(i,nps2mk)*((1.-SMALL)/(F(I,NPS1MK)+ | F(I,NPS2MK))) endif 41 CONTINUE RETURN END C