--- trunk/libf/dyn3d/ppm3d.f 2008/02/27 13:16:39 3 +++ trunk/libf/dyn3d/ppm3d.f 2010/04/06 17:52:58 32 @@ -80,10 +80,10 @@ C ============= C C Q(IMR,JNP,NLAY,NC): mixing ratios at current time (t) -C NC: total # of constituents -C IMR: first dimension (E-W); # of Grid intervals in E-W is IMR -C JNP: 2nd dimension (N-S); # of Grid intervals in N-S is JNP-1 -C NLAY: 3rd dimension (# of layers); vertical index increases from 1 at +C NC: total number of constituents +C IMR: first dimension (E-W); number of Grid intervals in E-W is IMR +C JNP: 2nd dimension (N-S); number of Grid intervals in N-S is JNP-1 +C NLAY: 3rd dimension (number of layers); vertical index increases from 1 at C the model top to NLAY near the surface (see fig. below). C It is assumed that 6 <= NLAY <= JNP (for dynamic memory allocation) C @@ -240,7 +240,7 @@ C winds are noisy near poles). C C Flux-Form Semi-Lagrangian transport in the East-West direction is used -C when and where Courant # is greater than one. +C when and where Courant number is greater than one. C C The user needs to change the parameter Jmax or Kmax if the resolution C is greater than 0.5 deg in N-S or 150 layers in the vertical direction. @@ -397,16 +397,13 @@ do J=2,JMR DTDX(j) = DT / ( DL*AE*COSP(J) ) -c print*,'dtdx=',dtdx(j) DTDX5(j) = 0.5*DTDX(j) enddo C DTDY = DT /(AE*DP) -c print*,'dtdy=',dtdy DTDY5 = 0.5*DTDY C -c write(6,*) 'J1=',J1,' J2=', J2 endif C C *********** End Initialization ********************** @@ -441,10 +438,10 @@ do 1500 k=1,NLAY C if(IGD.eq.0) then -C Convert winds on A-Grid to Courant # on C-Grid. +C Convert winds on A-Grid to Courant number on C-Grid. call A2C(U(1,1,k),V(1,1,k),IMR,JMR,j1,j2,CRX,CRY,dtdx5,DTDY5) else -C Convert winds on C-grid to Courant # +C Convert winds on C-grid to Courant number do 45 j=j1,j2 do 45 i=2,IMR 45 CRX(i,J) = dtdx(j)*U(i-1,j,k) @@ -726,7 +723,6 @@ DO j=1,JNP DO i=1,IMR Q(i,j,k,IC) = DQ(i,j,k,IC) / delp2(i,j,k) -c print*,'i=',i,'j=',j,'k=',k,'Q(i,j,k,IC)=',Q(i,j,k,IC) enddo enddo enddo @@ -863,18 +859,15 @@ AL(i,k) = wk1(i,k-1) + c1 + c2 * & ( wk2(i,k )*(c1*(A1 - A2)+A2*flux(i,k-1)) - & wk2(i,k-1)*A1*flux(i,k) ) -C print *,'AL1',i,k, AL(i,k) 12 CONTINUE 14 continue C do 20 i=1,IMR*NLAYM1 AR(i,1) = AL(i,2) -C print *,'AR1',i,AR(i,1) 20 continue C do 30 i=1,IMR*NLAY A6(i,1) = 3.*(wk1(i,1)+wk1(i,1) - (AL(i,1)+AR(i,1))) -C print *,'A61',i,A6(i,1) 30 continue C C****6***0*********0*********0*********0*********0*********0**********72 @@ -895,11 +888,8 @@ CM = wz2(i,1) / wk2(i,1) flux(i,2) = AR(i,1)+0.5*CM*(AL(i,1)-AR(i,1)+A6(i,1)*(1.-R23*CM)) else -C print *,'test2-0',i,j,wz2(i,1),wk2(i,2) CP= wz2(i,1) / wk2(i,2) -C print *,'testCP',CP flux(i,2) = AL(i,2)+0.5*CP*(AL(i,2)-AR(i,2)-A6(i,2)*(1.+R23*CP)) -C print *,'test2',i, AL(i,2),AR(i,2),A6(i,2),R23 endif 140 continue C @@ -1054,10 +1044,7 @@ c endif C LMT = IORD - 3 -c write(6,*) 'PPM option in E-W direction = ', LMT -c first = .false. -C endif -C + DO 10 i=1,IMR 10 AL(i) = 0.5*(p(i-1)+p(i)) + (DC(i-1) - DC(i))*R3 C @@ -1335,25 +1322,19 @@ wk(i,JNP+2) = p(i+IMH,JNP-2) wk(i+IMH,JNP+2) = p(i,JNP-2) enddo -c write(*,*) 'toto 1' -C -------------------------------- + IF(IAD.eq.2) then do j=j1-1,j2+1 do i=1,IMR -c write(*,*) 'avt NINT','i=',i,'j=',j JP = NINT(VA(i,j)) rv = JP - VA(i,j) -c write(*,*) 'VA=',VA(i,j), 'JP1=',JP,'rv=',rv JP = j - JP -c write(*,*) 'JP2=',JP a1 = 0.5*(wk(i,jp+1)+wk(i,jp-1)) - wk(i,jp) b1 = 0.5*(wk(i,jp+1)-wk(i,jp-1)) -c write(*,*) 'a1=',a1,'b1=',b1 ady(i,j) = wk(i,jp) + rv*(a1*rv + b1) - wk(i,j) enddo enddo -c write(*,*) 'toto 2' -C + ELSEIF(IAD.eq.1) then do j=j1-1,j2+1 do i=1,imr @@ -1980,22 +1961,3 @@ endif return end -C - subroutine zflip(q,im,km,nc) -C This routine flip the array q (in the vertical). - real q(im,km,nc) -C local dynamic array - real qtmp(im,km) -C - do 4000 IC = 1, nc -C - do 1000 k=1,km - do 1000 i=1,im - qtmp(i,k) = q(i,km+1-k,IC) -1000 continue -C - do 2000 i=1,im*km -2000 q(i,1,IC) = qtmp(i,1) -4000 continue - return - end