! ! $Header$ ! SUBROUTINE pentes_ini (q,w,masse,pbaru,pbarv,mode) USE comgeom2_mod_h USE comconst_mod, ONLY: pi, dtvr USE dimensions_mod, ONLY: iim, jjm, llm, ndm USE paramet_mod_h IMPLICIT NONE !======================================================================= ! Adaptation LMDZ: A.Armengaud (LGGE) ! ---------------- ! ! ******************************************************************** ! Transport des traceurs par la methode des pentes ! ******************************************************************** ! Reference possible : Russel. G.L., Lerner J.A.: ! A new Finite-Differencing Scheme for Traceur Transport ! Equation , Journal of Applied Meteorology, pp 1483-1498,dec. 81 ! ******************************************************************** ! q,w,masse,pbaru et pbarv ! sont des arguments d'entree pour le s-pg .... ! !======================================================================= ! Arguments: ! ---------- integer :: mode REAL :: pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm ) REAL :: q( iip1,jjp1,llm,0:3) REAL :: w( ip1jmp1,llm ) REAL :: masse( iip1,jjp1,llm) ! Local: ! ------ LOGICAL :: limit REAL :: sm ( iip1,jjp1, llm ) REAL :: s0( iip1,jjp1,llm ), sx( iip1,jjp1,llm ) REAL :: sy( iip1,jjp1,llm ), sz( iip1,jjp1,llm ) real :: masn,mass,zz INTEGER :: i,j,l,iq ! modif Fred 24 03 96 real :: sinlon(iip1),sinlondlon(iip1) real :: coslon(iip1),coslondlon(iip1) save sinlon,coslon,sinlondlon,coslondlon real :: dyn1,dyn2,dys1,dys2 real :: qpn,qps,dqzpn,dqzps real :: smn,sms,s0n,s0s,sxn(iip1),sxs(iip1) real :: qmin,zq,pente_max ! REAL :: SSUM integer :: ismax,ismin,lati,latf EXTERNAL SSUM, ismin,ismax logical :: first save first ! fin modif ! EXTERNAL masskg EXTERNAL advx EXTERNAL advy EXTERNAL advz ! modif Fred 24 03 96 data first/.true./ limit = .TRUE. pente_max=2 ! if (mode.eq.1.or.mode.eq.3) then ! if (mode.eq.1) then if (mode.ge.1) then lati=2 latf=jjm else lati=1 latf=jjp1 endif qmin=0.4995 qmin=0. if(first) then print*,'SCHEMA AMONT NOUVEAU' first=.false. do i=2,iip1 coslon(i)=cos(rlonv(i)) sinlon(i)=sin(rlonv(i)) coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi print*,coslondlon(i),sinlondlon(i) enddo coslon(1)=coslon(iip1) coslondlon(1)=coslondlon(iip1) sinlon(1)=sinlon(iip1) sinlondlon(1)=sinlondlon(iip1) print*,'sum sinlondlon ',ssum(iim,sinlondlon,1)/sinlondlon(1) print*,'sum coslondlon ',ssum(iim,coslondlon,1)/coslondlon(1) DO l = 1,llm DO j = 1,jjp1 DO i = 1,iip1 q ( i,j,l,1 )=0. q ( i,j,l,2 )=0. q ( i,j,l,3 )=0. ENDDO ENDDO ENDDO endif ! Fin modif Fred ! *** q contient les qqtes de traceur avant l'advection ! *** Affectation des tableaux S a partir de Q ! *** Rem : utilisation de SCOPY ulterieurement DO l = 1,llm DO j = 1,jjp1 DO i = 1,iip1 s0( i,j,llm+1-l ) = q ( i,j,l,0 ) sx( i,j,llm+1-l ) = q ( i,j,l,1 ) sy( i,j,llm+1-l ) = q ( i,j,l,2 ) sz( i,j,llm+1-l ) = q ( i,j,l,3 ) ENDDO ENDDO ENDDO ! PRINT*,'----- S0 just before conversion -------' ! PRINT*,'S0(16,12,1)=',s0(16,12,1) ! PRINT*,'Q(16,12,1,4)=',q(16,12,1,4) ! *** On calcule la masse d'air en kg DO l = 1,llm DO j = 1,jjp1 DO i = 1,iip1 sm ( i,j,llm+1-l)=masse( i,j,l ) ENDDO ENDDO ENDDO ! *** On converti les champs S en atome (resp. kg) ! *** Les routines d'advection traitent les champs ! *** a advecter si ces derniers sont en atome (resp. kg) ! *** A optimiser !!! DO l = 1,llm DO j = 1,jjp1 DO i = 1,iip1 s0(i,j,l) = s0(i,j,l) * sm ( i,j,l ) sx(i,j,l) = sx(i,j,l) * sm ( i,j,l ) sy(i,j,l) = sy(i,j,l) * sm ( i,j,l ) sz(i,j,l) = sz(i,j,l) * sm ( i,j,l ) ENDDO ENDDO ENDDO ! ss0 = 0. ! DO l = 1,llm ! DO j = 1,jjp1 ! DO i = 1,iim ! ss0 = ss0 + s0 ( i,j,l ) ! ENDDO ! ENDDO ! ENDDO ! PRINT*, 'valeur tot s0 avant advection=',ss0 ! *** Appel des subroutines d'advection en X, en Y et en Z ! *** Advection avec "time-splitting" !----------------------------------------------------------- ! PRINT*,'----- S0 just before ADVX -------' ! PRINT*,'S0(16,12,1)=',s0(16,12,1) !----------------------------------------------------------- ! do l=1,llm ! do j=1,jjp1 ! do i=1,iip1 ! zq=s0(i,j,l)/sm(i,j,l) ! if(zq.lt.qmin) ! , print*,'avant advx1, s0(',i,',',j,',',l,')=',zq ! enddo ! enddo ! enddo !CC if(mode.eq.2) then do l=1,llm s0s=0. s0n=0. dyn1=0. dys1=0. dyn2=0. dys2=0. smn=0. sms=0. do i=1,iim smn=smn+sm(i,1,l) sms=sms+sm(i,jjp1,l) s0n=s0n+s0(i,1,l) s0s=s0s+s0(i,jjp1,l) zz=sy(i,1,l)/sm(i,1,l) dyn1=dyn1+sinlondlon(i)*zz dyn2=dyn2+coslondlon(i)*zz zz=sy(i,jjp1,l)/sm(i,jjp1,l) dys1=dys1+sinlondlon(i)*zz dys2=dys2+coslondlon(i)*zz enddo do i=1,iim sy(i,1,l)=dyn1*sinlon(i)+dyn2*coslon(i) sy(i,jjp1,l)=dys1*sinlon(i)+dys2*coslon(i) enddo do i=1,iim s0(i,1,l)=s0n/smn+sy(i,1,l) s0(i,jjp1,l)=s0s/sms-sy(i,jjp1,l) enddo s0(iip1,1,l)=s0(1,1,l) s0(iip1,jjp1,l)=s0(1,jjp1,l) do i=1,iim sxn(i)=s0(i+1,1,l)-s0(i,1,l) sxs(i)=s0(i+1,jjp1,l)-s0(i,jjp1,l) ! on rerentre les masses enddo do i=1,iim sy(i,1,l)=sy(i,1,l)*sm(i,1,l) sy(i,jjp1,l)=sy(i,jjp1,l)*sm(i,jjp1,l) s0(i,1,l)=s0(i,1,l)*sm(i,1,l) s0(i,jjp1,l)=s0(i,jjp1,l)*sm(i,jjp1,l) enddo sxn(iip1)=sxn(1) sxs(iip1)=sxs(1) do i=1,iim sx(i+1,1,l)=0.25*(sxn(i)+sxn(i+1))*sm(i+1,1,l) sx(i+1,jjp1,l)=0.25*(sxs(i)+sxs(i+1))*sm(i+1,jjp1,l) enddo s0(iip1,1,l)=s0(1,1,l) s0(iip1,jjp1,l)=s0(1,jjp1,l) sy(iip1,1,l)=sy(1,1,l) sy(iip1,jjp1,l)=sy(1,jjp1,l) sx(1,1,l)=sx(iip1,1,l) sx(1,jjp1,l)=sx(iip1,jjp1,l) enddo endif if (mode.eq.4) then do l=1,llm do i=1,iip1 sx(i,1,l)=0. sx(i,jjp1,l)=0. sy(i,1,l)=0. sy(i,jjp1,l)=0. enddo enddo endif call limx(s0,sx,sm,pente_max) ! call minmaxq(zq,1.e33,-1.e33,'avant advx ') call advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf) ! call minmaxq(zq,1.e33,-1.e33,'avant advy ') if (mode.eq.4) then do l=1,llm do i=1,iip1 sx(i,1,l)=0. sx(i,jjp1,l)=0. sy(i,1,l)=0. sy(i,jjp1,l)=0. enddo enddo endif call limy(s0,sy,sm,pente_max) call advy( limit,.5*dtvr,pbarv,sm,s0,sx,sy,sz ) ! call minmaxq(zq,1.e33,-1.e33,'avant advz ') do j=1,jjp1 do i=1,iip1 sz(i,j,1)=0. sz(i,j,llm)=0. enddo enddo call limz(s0,sz,sm,pente_max) call advz( limit,dtvr,w,sm,s0,sx,sy,sz ) if (mode.eq.4) then do l=1,llm do i=1,iip1 sx(i,1,l)=0. sx(i,jjp1,l)=0. sy(i,1,l)=0. sy(i,jjp1,l)=0. enddo enddo endif call limy(s0,sy,sm,pente_max) call advy( limit,.5*dtvr,pbarv,sm,s0,sx,sy,sz ) do l=1,llm do j=1,jjp1 sm(iip1,j,l)=sm(1,j,l) s0(iip1,j,l)=s0(1,j,l) sx(iip1,j,l)=sx(1,j,l) sy(iip1,j,l)=sy(1,j,l) sz(iip1,j,l)=sz(1,j,l) enddo enddo ! call minmaxq(zq,1.e33,-1.e33,'avant advx ') if (mode.eq.4) then do l=1,llm do i=1,iip1 sx(i,1,l)=0. sx(i,jjp1,l)=0. sy(i,1,l)=0. sy(i,jjp1,l)=0. enddo enddo endif call limx(s0,sx,sm,pente_max) call advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf) ! call minmaxq(zq,1.e33,-1.e33,'apres advx ') ! do l=1,llm ! do j=1,jjp1 ! do i=1,iip1 ! zq=s0(i,j,l)/sm(i,j,l) ! if(zq.lt.qmin) ! , print*,'apres advx2, s0(',i,',',j,',',l,')=',zq ! enddo ! enddo ! enddo ! *** On repasse les S dans la variable q directement 14/10/94 ! On revient a des rapports de melange en divisant par la masse ! En dehors des poles: DO l = 1,llm DO j = 1,jjp1 DO i = 1,iim q(i,j,llm+1-l,0)=s0(i,j,l)/sm(i,j,l) q(i,j,llm+1-l,1)=sx(i,j,l)/sm(i,j,l) q(i,j,llm+1-l,2)=sy(i,j,l)/sm(i,j,l) q(i,j,llm+1-l,3)=sz(i,j,l)/sm(i,j,l) ENDDO ENDDO ENDDO ! Traitements specifiques au pole if(mode.ge.1) then DO l=1,llm ! filtrages aux poles masn=ssum(iim,sm(1,1,l),1) mass=ssum(iim,sm(1,jjp1,l),1) qpn=ssum(iim,s0(1,1,l),1)/masn qps=ssum(iim,s0(1,jjp1,l),1)/mass dqzpn=ssum(iim,sz(1,1,l),1)/masn dqzps=ssum(iim,sz(1,jjp1,l),1)/mass do i=1,iip1 q( i,1,llm+1-l,3)=dqzpn q( i,jjp1,llm+1-l,3)=dqzps q( i,1,llm+1-l,0)=qpn q( i,jjp1,llm+1-l,0)=qps enddo if(mode.eq.3) then dyn1=0. dys1=0. dyn2=0. dys2=0. do i=1,iim dyn1=dyn1+sinlondlon(i)*sy(i,1,l)/sm(i,1,l) dyn2=dyn2+coslondlon(i)*sy(i,1,l)/sm(i,1,l) dys1=dys1+sinlondlon(i)*sy(i,jjp1,l)/sm(i,jjp1,l) dys2=dys2+coslondlon(i)*sy(i,jjp1,l)/sm(i,jjp1,l) enddo do i=1,iim q(i,1,llm+1-l,2)= & (sinlon(i)*dyn1+coslon(i)*dyn2) q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)+q(i,1,llm+1-l,2) q(i,jjp1,llm+1-l,2)= & (sinlon(i)*dys1+coslon(i)*dys2) q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0) & -q(i,jjp1,llm+1-l,2) enddo endif if(mode.eq.1) then ! on filtre les valeurs au bord de la "grande maille pole" dyn1=0. dys1=0. dyn2=0. dys2=0. do i=1,iim zz=s0(i,2,l)/sm(i,2,l)-q(i,1,llm+1-l,0) dyn1=dyn1+sinlondlon(i)*zz dyn2=dyn2+coslondlon(i)*zz zz=q(i,jjp1,llm+1-l,0)-s0(i,jjm,l)/sm(i,jjm,l) dys1=dys1+sinlondlon(i)*zz dys2=dys2+coslondlon(i)*zz enddo do i=1,iim q(i,1,llm+1-l,2)= & (sinlon(i)*dyn1+coslon(i)*dyn2)/2. q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)+q(i,1,llm+1-l,2) q(i,jjp1,llm+1-l,2)= & (sinlon(i)*dys1+coslon(i)*dys2)/2. q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0) & -q(i,jjp1,llm+1-l,2) enddo q(iip1,1,llm+1-l,0)=q(1,1,llm+1-l,0) q(iip1,jjp1,llm+1-l,0)=q(1,jjp1,llm+1-l,0) do i=1,iim sxn(i)=q(i+1,1,llm+1-l,0)-q(i,1,llm+1-l,0) sxs(i)=q(i+1,jjp1,llm+1-l,0)-q(i,jjp1,llm+1-l,0) enddo sxn(iip1)=sxn(1) sxs(iip1)=sxs(1) do i=1,iim q(i+1,1,llm+1-l,1)=0.25*(sxn(i)+sxn(i+1)) q(i+1,jjp1,llm+1-l,1)=0.25*(sxs(i)+sxs(i+1)) enddo q(1,1,llm+1-l,1)=q(iip1,1,llm+1-l,1) q(1,jjp1,llm+1-l,1)=q(iip1,jjp1,llm+1-l,1) endif ENDDO endif ! bouclage en longitude do iq=0,3 do l=1,llm do j=1,jjp1 q(iip1,j,l,iq)=q(1,j,l,iq) enddo enddo enddo ! PRINT*, ' SORTIE DE PENTES --- ca peut glisser ....' DO l = 1,llm DO j = 1,jjp1 DO i = 1,iip1 IF (q(i,j,l,0).lt.0.) THEN ! PRINT*,'------------ BIP-----------' ! PRINT*,'Q0(',i,j,l,')=',q(i,j,l,0) ! PRINT*,'QX(',i,j,l,')=',q(i,j,l,1) ! PRINT*,'QY(',i,j,l,')=',q(i,j,l,2) ! PRINT*,'QZ(',i,j,l,')=',q(i,j,l,3) ! PRINT*,' PBL EN SORTIE DE PENTES' q(i,j,l,0)=0. ! STOP ENDIF ENDDO ENDDO ENDDO ! PRINT*, '-------------------------------------------' do l=1,llm do j=1,jjp1 do i=1,iip1 if(q(i,j,l,0).lt.qmin) & print*,'apres pentes, s0(',i,',',j,',',l,')=',q(i,j,l,0) enddo enddo enddo RETURN END SUBROUTINE pentes_ini