Index: /LMDZ.3.3/branches/rel-LF/libf/dyn3d/gcm.F =================================================================== --- /LMDZ.3.3/branches/rel-LF/libf/dyn3d/gcm.F (revision 494) +++ /LMDZ.3.3/branches/rel-LF/libf/dyn3d/gcm.F (revision 495) @@ -153,4 +153,5 @@ LOGICAL ok_nudge PARAMETER (ok_nudge = .false.) +c c----------------------------------------------------------------------- @@ -367,4 +368,5 @@ c ---------------------------------- + 1 CONTINUE @@ -449,12 +451,17 @@ ELSE IF( iq.EQ. nqmx ) THEN c - iapp_tracvl = 5 +c iapp_tracvl = 5 + iapp_tracvl = iperiod + print*,'***WARNING***: iapp_tracvl = iperiod' c cccc iapp_tracvl est la frequence en pas du groupement des flux cccc de masse pour Van-Leer dans la routine tracvl . c - CALL vanleer(numvanle,iapp_tracvl,nqmx,q,pbaru,pbarv, +c CALL vanleer(numvanle,iapp_tracvl,nqmx,q,pbaru,pbarv, + CALL vanleer(numvanle,iapp_tracvl,2,q,pbaru,pbarv, * p, masse, dq, iadv(1), teta, pk ) c + print*,'***WARNING***: Appel vanleer avec 2 traceurs' + print*,' et non nqmx' c ... Modif F.Codron .... c @@ -483,4 +490,5 @@ $ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis , $ finvmaold ) + c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd) @@ -660,4 +668,5 @@ IF( itau. EQ. itaufinp1 ) then + abort_message = 'Simulation finished' call abort_gcm(modname,abort_message,0) Index: /LMDZ.3.3/branches/rel-LF/libf/filtrez/filtreg.F =================================================================== --- /LMDZ.3.3/branches/rel-LF/libf/filtrez/filtreg.F (revision 494) +++ /LMDZ.3.3/branches/rel-LF/libf/filtrez/filtreg.F (revision 495) @@ -49,5 +49,5 @@ INTEGER i,j,l,k INTEGER iim2,immjm - INTEGER jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil + INTEGER jdfil1,jdfil2,jffil1,jffil2,jdfil(2),jffil(2) REAL champ( iip1,nlat,nbniv) @@ -56,5 +56,6 @@ , , matricevn(iim,iim,nfilvn),matricevs(iim,iim,nfilvs) , , matrinvn(iim,iim,nfilun),matrinvs (iim,iim,nfilus) - REAL eignq(iim), sdd1(iim),sdd2(iim) +cIM REAL eignq(iim), sdd1(iim),sdd2(iim) + REAL eignq(iim,nlat,nbniv), sdd1(iim),sdd2(iim) LOGICAL griscal INTEGER hemisph, iaire @@ -127,18 +128,22 @@ END IF c + jdfil(1) = jdfil1 + jffil(1) = jffil1 + jdfil(2) = jdfil2 + jffil(2) = jffil2 c DO 100 hemisph = 1, 2 c - IF ( hemisph.EQ.1 ) THEN - jdfil = jdfil1 - jffil = jffil1 - ELSE - jdfil = jdfil2 - jffil = jffil2 - END IF +c IF ( hemisph.EQ.1 ) THEN +c jdfil = jdfil1 +c jffil = jffil1 +c ELSE +c jdfil = jdfil2 +c jffil = jffil2 +c END IF DO 50 l = 1, nbniv - DO 30 j = jdfil,jffil + DO 30 j = jdfil(hemisph),jffil(hemisph) @@ -147,4 +152,6 @@ 5 CONTINUE c + 30 CONTINUE + 50 CONTINUE IF( hemisph. EQ. 1 ) THEN @@ -152,57 +159,99 @@ IF( ifiltre. EQ. -2 ) THEN #ifdef CRAY + DO l = 1, nbniv + DO j = jdfil(hemisph),jffil(hemisph) CALL MXVA( matrinvn(1,1,j), 1, iim, champ(1,j,l), 1, eignq , * 1, iim, iim ) -#else -#ifdef BLAS + ENDDO + ENDDO +#else +#ifdef BLAS + DO l = 1, nbniv + DO j = jdfil(hemisph),jffil(hemisph) CALL SGEMV("N", iim,iim, 1.0, matrinvn(1,1,j),iim, . champ(1,j,l), 1, 0.0, eignq, 1) -#else - DO k = 1, iim - eignq(k) = 0.0 - ENDDO - DO k = 1, iim - DO i = 1, iim - eignq(k) = eignq(k) + matrinvn(k,i,j)*champ(i,j,l) - ENDDO - ENDDO + ENDDO + ENDDO +#else + DO l = 1, nbniv + DO j = jdfil(hemisph),jffil(hemisph) + DO k = 1, iim +c eignq(k) = 0.0 + eignq(k,j,l) = 0.0 + ENDDO + DO k = 1, iim + DO i = 1, iim +c eignq(k) = eignq(k) + matrinvn(k,i,j)*champ(i,j,l) + eignq(k,j,l) = eignq(k,j,l) + matrinvn(k,i,j)*champ(i,j,l) + ENDDO + ENDDO + ENDDO + ENDDO #endif #endif ELSE IF ( griscal ) THEN #ifdef CRAY + DO l = 1, nbniv + DO j = jdfil(hemisph),jffil(hemisph) CALL MXVA( matriceun(1,1,j), 1, iim, champ(1,j,l), 1, eignq , * 1, iim, iim ) -#else -#ifdef BLAS + ENDDO + ENDDO +#else +#ifdef BLAS + DO l = 1, nbniv + DO j = jdfil(hemisph),jffil(hemisph) CALL SGEMV("N", iim,iim, 1.0, matriceun(1,1,j),iim, . champ(1,j,l), 1, 0.0, eignq, 1) -#else - DO k = 1, iim - eignq(k) = 0.0 - ENDDO - DO i = 1, iim - DO k = 1, iim - eignq(k) = eignq(k) + matriceun(k,i,j)*champ(i,j,l) - ENDDO - ENDDO + ENDDO + ENDDO +#else + DO l = 1, nbniv + DO j = jdfil(hemisph),jffil(hemisph) + DO k = 1, iim +c eignq(k) = 0.0 + eignq(k,j,l) = 0.0 + ENDDO + DO i = 1, iim + DO k = 1, iim +c eignq(k) = eignq(k) + matriceun(k,i,j)*champ(i,j,l) + eignq(k,j,l) = eignq(k,j,l) + matriceun(k,i,j)*champ(i,j,l) + ENDDO + ENDDO + ENDDO + ENDDO #endif #endif ELSE #ifdef CRAY + DO l = 1, nbniv + DO j = jdfil(hemisph),jffil(hemisph) CALL MXVA( matricevn(1,1,j), 1, iim, champ(1,j,l), 1, eignq , * 1, iim, iim ) -#else -#ifdef BLAS + ENDDO + ENDDO +#else +#ifdef BLAS + DO l = 1, nbniv + DO j = jdfil(hemisph),jffil(hemisph) CALL SGEMV("N", iim,iim, 1.0, matricevn(1,1,j),iim, . champ(1,j,l), 1, 0.0, eignq, 1) -#else - DO k = 1, iim - eignq(k) = 0.0 - ENDDO - DO i = 1, iim - DO k = 1, iim - eignq(k) = eignq(k) + matricevn(k,i,j)*champ(i,j,l) - ENDDO - ENDDO + ENDDO + ENDDO +#else + DO l = 1, nbniv + DO j = jdfil(hemisph),jffil(hemisph) + DO k = 1, iim +c eignq(k) = 0.0 + eignq(k,j,l) = 0.0 + ENDDO + DO i = 1, iim + DO k = 1, iim +c eignq(k) = eignq(k) + matricevn(k,i,j)*champ(i,j,l) + eignq(k,j,l) = eignq(k,j,l) + matricevn(k,i,j)*champ(i,j,l) + ENDDO + ENDDO + ENDDO + ENDDO #endif #endif @@ -213,57 +262,102 @@ IF( ifiltre. EQ. -2 ) THEN #ifdef CRAY + DO l = 1, nbniv + DO j = jdfil(hemisph),jffil(hemisph) CALL MXVA( matrinvs(1,1,j-jfiltsu+1), 1, iim, champ(1,j,l),1 , * eignq, 1, iim, iim ) -#else -#ifdef BLAS + ENDDO + ENDDO +#else +#ifdef BLAS + DO l = 1, nbniv + DO j = jdfil(hemisph),jffil(hemisph) CALL SGEMV("N", iim,iim, 1.0, matrinvs(1,1,j-jfiltsu+1),iim, . champ(1,j,l), 1, 0.0, eignq, 1) -#else - DO k = 1, iim - eignq(k) = 0.0 - ENDDO - DO i = 1, iim - DO k = 1, iim - eignq(k) = eignq(k) + matrinvs(k,i,j-jfiltsu+1)*champ(i,j,l) - ENDDO - ENDDO + ENDDO + ENDDO +#else + DO l = 1, nbniv + DO j = jdfil(hemisph),jffil(hemisph) + DO k = 1, iim +c eignq(k) = 0.0 + eignq(k,j,l) = 0.0 + ENDDO + DO i = 1, iim + DO k = 1, iim +c eignq(k) = eignq(k) + matrinvs(k,i,j-jfiltsu+1)*champ(i,j,l) + eignq(k,j,l) = eignq(k,j,l) + + .matrinvs(k,i,j-jfiltsu+1)*champ(i,j,l) + ENDDO + ENDDO + ENDDO + ENDDO #endif #endif ELSE IF ( griscal ) THEN #ifdef CRAY + DO l = 1, nbniv + DO j = jdfil(hemisph),jffil(hemisph) CALL MXVA( matriceus(1,1,j-jfiltsu+1), 1, iim, champ(1,j,l),1 , * eignq, 1, iim, iim ) -#else -#ifdef BLAS + ENDDO + ENDDO +#else +#ifdef BLAS + DO l = 1, nbniv + DO j = jdfil(hemisph),jffil(hemisph) CALL SGEMV("N", iim,iim, 1.0, matriceus(1,1,j-jfiltsu+1),iim, . champ(1,j,l), 1, 0.0, eignq, 1) -#else - DO k = 1, iim - eignq(k) = 0.0 - ENDDO - DO i = 1, iim - DO k = 1, iim - eignq(k) = eignq(k) + matriceus(k,i,j-jfiltsu+1)*champ(i,j,l) - ENDDO - ENDDO + ENDDO + ENDDO +#else + DO l = 1, nbniv + DO j = jdfil(hemisph),jffil(hemisph) + DO k = 1, iim +c eignq(k) = 0.0 + eignq(k,j,l) = 0.0 + ENDDO + DO i = 1, iim + DO k = 1, iim +c eignq(k) = eignq(k) + matriceus(k,i,j-jfiltsu+1)*champ(i,j,l) + eignq(k,j,l) = eignq(k,j,l) + + .matriceus(k,i,j-jfiltsu+1)*champ(i,j,l) + ENDDO + ENDDO + ENDDO + ENDDO #endif #endif ELSE #ifdef CRAY + DO l = 1, nbniv + DO j = jdfil(hemisph),jffil(hemisph) CALL MXVA( matricevs(1,1,j-jfiltsv+1), 1, iim, champ(1,j,l),1 , * eignq, 1, iim, iim ) -#else -#ifdef BLAS + ENDDO + ENDDO +#else +#ifdef BLAS + DO l = 1, nbniv + DO j = jdfil(hemisph),jffil(hemisph) CALL SGEMV("N", iim,iim, 1.0, matricevs(1,1,j-jfiltsv+1),iim, . champ(1,j,l), 1, 0.0, eignq, 1) -#else - DO k = 1, iim - eignq(k) = 0.0 - ENDDO - DO i = 1, iim - DO k = 1, iim - eignq(k) = eignq(k) + matricevs(k,i,j-jfiltsv+1)*champ(i,j,l) - ENDDO - ENDDO + ENDDO + ENDDO +#else + DO l = 1, nbniv + DO j = jdfil(hemisph),jffil(hemisph) + DO k = 1, iim +c eignq(k) = 0.0 + eignq(k,j,l) = 0.0 + ENDDO + DO i = 1, iim + DO k = 1, iim +c eignq(k) = eignq(k) + matricevs(k,i,j-jfiltsv+1)*champ(i,j,l) + eignq(k,j,l) = eignq(k,j,l) + + .matricevs(k,i,j-jfiltsv+1)*champ(i,j,l) + ENDDO + ENDDO + ENDDO + ENDDO #endif #endif @@ -273,18 +367,32 @@ c IF( ifiltre.EQ. 2 ) THEN + DO l = 1, nbniv + DO j = jdfil(hemisph),jffil(hemisph) DO 15 i = 1, iim - champ( i,j,l ) = ( champ(i,j,l) + eignq(i) ) * sdd2(i) +c champ( i,j,l ) = ( champ(i,j,l) + eignq(i) ) * sdd2(i) + champ( i,j,l ) = ( champ(i,j,l) + eignq(i,j,l) ) * sdd2(i) 15 CONTINUE + ENDDO + ENDDO ELSE + DO l = 1, nbniv + DO j = jdfil(hemisph),jffil(hemisph) DO 16 i=1,iim - champ( i,j,l ) = ( champ(i,j,l) - eignq(i) ) * sdd2(i) +c champ( i,j,l ) = ( champ(i,j,l) - eignq(i) ) * sdd2(i) + champ( i,j,l ) = ( champ(i,j,l) - eignq(i,j,l) ) * sdd2(i) 16 CONTINUE + ENDDO + ENDDO ENDIF c + DO l = 1, nbniv + DO j = jdfil(hemisph),jffil(hemisph) champ( iip1,j,l ) = champ( 1,j,l ) -c - 30 CONTINUE -c - 50 CONTINUE + ENDDO + ENDDO +c +c 30 CONTINUE +c +c 50 CONTINUE c 100 CONTINUE Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/clouds_gno.F =================================================================== --- /LMDZ.3.3/branches/rel-LF/libf/phylmd/clouds_gno.F (revision 494) +++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/clouds_gno.F (revision 495) @@ -48,7 +48,8 @@ REAL xx(klon), aux(klon), coeff(klon), block(klon) REAL dist(klon), fprime(klon), det(klon) - REAL pi, u(klon), v(klon), erfu(klon), erfv(klon) + REAL pi, u(klon), v(klon), erfcu(klon), erfcv(klon) REAL xx1(klon), xx2(klon) real erf,kkk + real sqrtpi,sqrt2,zx1,zx2,exdel c lconv = true si le calcul a converge (entre autre si qsub < min_q) LOGICAL lconv(klon) @@ -56,4 +57,6 @@ pi = ACOS(-1.) + sqrtpi=sqrt(pi) + sqrt2=sqrt(2.) ptconv=.false. @@ -123,6 +126,8 @@ xx(i) = -0.0001 else - xx1(i)=-SQRT(2.)*vmax(i)*(1.0-SQRT(1.0+delta(i)/(vmax(i)**2.))) - xx2(i)=-SQRT(2.)*vmax(i)*(1.0+SQRT(1.0+delta(i)/(vmax(i)**2.))) + zx1=-sqrt2*vmax(i) + zx2=SQRT(1.0+delta(i)/(vmax(i)**2.)) + xx1(i)=zx1*(1.0-zx2) + xx2(i)=zx1*(1.0+zx2) xx(i) = 1.01 * xx1(i) if ( xx1(i) .GE. 0.0 ) xx(i) = 0.5*xx2(i) @@ -143,6 +148,6 @@ if (.not.lconv(i)) then - u(i) = delta(i)/(xx(i)*sqrt(2.)) + xx(i)/(2.*sqrt(2.)) - v(i) = delta(i)/(xx(i)*sqrt(2.)) - xx(i)/(2.*sqrt(2.)) + u(i) = delta(i)/(xx(i)*sqrt2) + xx(i)/(2.*sqrt2) + v(i) = delta(i)/(xx(i)*sqrt2) - xx(i)/(2.*sqrt2) IF ( v(i) .GT. vmax(i) ) THEN @@ -153,13 +158,13 @@ c -- use asymptotic expression of erf for u and v large: c ( -> analytic solution for xx ) - - aux(i) = 2.0*delta(i)*(1.-beta(i)*EXP(delta(i))) - : /(1.+beta(i)*EXP(delta(i))) + exdel=beta(i)*EXP(delta(i)) + aux(i) = 2.0*delta(i)*(1.-exdel) + : /(1.+exdel) if (aux(i).lt.0.) then - print*,'AUX(',i,',',k,')<0',aux(i),delta(i),beta(i) +c print*,'AUX(',i,',',k,')<0',aux(i),delta(i),beta(i) aux(i)=0. endif xx(i) = -SQRT(aux(i)) - block(i) = EXP(-v(i)*v(i)) / v(i) / SQRT(pi) + block(i) = EXP(-v(i)*v(i)) / v(i) / sqrtpi dist(i) = 0.0 fprime(i) = 1.0 @@ -169,9 +174,9 @@ c -- erfv -> 1.0, use an asymptotic expression of erfv for v large: - erfu(i) = ERF(u(i)) + erfcu(i) = 1.0-ERF(u(i)) c !!! ATTENTION : rajout d'un seuil pour l'exponentiel - aux(i) = SQRT(pi)*(1.0-erfu(i))*EXP(min(v(i)*v(i),100.)) + aux(i) = sqrtpi*erfcu(i)*EXP(min(v(i)*v(i),100.)) coeff(i) = 1.0 - 1./2./(v(i)**2.) + 3./4./(v(i)**4.) - block(i) = coeff(i) * EXP(-v(i)*v(i)) / v(i) / SQRT(pi) + block(i) = coeff(i) * EXP(-v(i)*v(i)) / v(i) / sqrtpi dist(i) = v(i) * aux(i) / coeff(i) - beta(i) fprime(i) = 2.0 / xx(i) * (v(i)**2.) @@ -185,23 +190,23 @@ c -- general case: - erfu(i) = ERF(u(i)) - erfv(i) = ERF(v(i)) - block(i) = 1.0-erfv(i) - dist(i) = (1.0 - erfu(i)) / (1.0 - erfv(i)) - beta(i) + erfcu(i) = 1.0-ERF(u(i)) + erfcv(i) = 1.0-ERF(v(i)) + block(i) = erfcv(i) + dist(i) = erfcu(i) / erfcv(i) - beta(i) zu2(i)=u(i)*u(i) zv2(i)=v(i)*v(i) if(zu2(i).gt.20..or. zv2(i).gt.20.) then - print*,'ATTENTION !!! xx(',i,') =', xx(i) - print*,'ATTENTION !!! klon,ND,R,RS,QSUB,PTCONV,RATQSC,CLDF', - .klon,ND,R(i,k),RS(i,k),QSUB(i,k),PTCONV(i,k),RATQSC(i,k), - .CLDF(i,k) - print*,'ATTENTION !!! zu2 zv2 =',zu2(i),zv2(i) +c print*,'ATTENTION !!! xx(',i,') =', xx(i) +c print*,'ATTENTION !!! klon,ND,R,RS,QSUB,PTCONV,RATQSC,CLDF', +c .klon,ND,R(i,k),RS(i,k),QSUB(i,k),PTCONV(i,k),RATQSC(i,k), +c .CLDF(i,k) +c print*,'ATTENTION !!! zu2 zv2 =',zu2(i),zv2(i) zu2(i)=20. zv2(i)=20. fprime(i) = 0. else - fprime(i) = 2. /SQRT(pi) /xx(i) /(1.0-erfv(i))**2. - : * ( (1.0-erfv(i))*v(i)*EXP(-zu2(i)) - : - (1.0-erfu(i))*u(i)*EXP(-zv2(i)) ) + fprime(i) = 2. /sqrtpi /xx(i) /erfcv(i)**2. + : * ( erfcv(i)*v(i)*EXP(-zu2(i)) + : - erfcu(i)*u(i)*EXP(-zv2(i)) ) endif ENDIF ! x @@ -219,14 +224,4 @@ ratqsc(i,k)=sqrt(exp(ratqsc(i,k))-1.) CLDF(i,K) = 0.5 * block(i) -cccccccccccccccccccccccc -c kkk=-sqrt(log(1.+ratqsc(i,k)**2)) -c u(i)=delta(i)/(kkk*sqrt(2.))-kkk/(2.*sqrt(2.)) -c v(i)=delta(i)/(kkk*sqrt(2.))+kkk/(2.*sqrt(2.)) -c erfu(i)=erf(u(i)) -c erfv(i)=erf(v(i)) -c print*,'SIG ',k,qsub(i,k) -c s ,mu(i)*((1.-erfv(i))/(1.-erfu(i))-qsat(i)/mu(i)) -c s ,0.5*erfu(i) -cccccccccccccccccccccccc else xx(i) = xx(i) - dist(i)/fprime(i) Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/cv3_routines.F =================================================================== --- /LMDZ.3.3/branches/rel-LF/libf/phylmd/cv3_routines.F (revision 494) +++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/cv3_routines.F (revision 495) @@ -804,16 +804,15 @@ 110 continue - do 121 j=1,ntra -ccccc do 111 k=1,nl+1 - do 111 k=1,nd - nn=0 - do 101 i=1,len - if(iflag1(i).eq.0)then - nn=nn+1 - tra(nn,k,j)=tra1(i,k,j) - endif - 101 continue - 111 continue - 121 continue +c do 121 j=1,ntra +c do 111 k=1,nd +c nn=0 +c do 101 i=1,len +c if(iflag1(i).eq.0)then +c nn=nn+1 +c tra(nn,k,j)=tra1(i,k,j) +c endif +c 101 continue +c 111 continue +c 121 continue if (nn.ne.ncum) then @@ -1493,13 +1492,13 @@ 400 continue - do k=1,ntra - do j=1,nd ! instead nlp - do i=1,nd ! instead nlp - do il=1,ncum - traent(il,i,j,k)=tra(il,j,k) - enddo - enddo - enddo - enddo +c do k=1,ntra +c do j=1,nd ! instead nlp +c do i=1,nd ! instead nlp +c do il=1,ncum +c traent(il,i,j,k)=tra(il,j,k) +c enddo +c enddo +c enddo +c enddo zm(:,:)=0. @@ -1557,15 +1556,15 @@ 710 continue - do k=1,ntra - do j=minorig,nl - do il=1,ncum - if( (i.ge.icb(il)).and.(i.le.inb(il)).and. - : (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then - traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k) - : +(1.-sij(il,i,j))*tra(il,nk(il),k) - endif - enddo - enddo - enddo +c do k=1,ntra +c do j=minorig,nl +c do il=1,ncum +c if( (i.ge.icb(il)).and.(i.le.inb(il)).and. +c : (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then +c traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k) +c : +(1.-sij(il,i,j))*tra(il,nk(il),k) +c endif +c enddo +c enddo +c enddo c @@ -1590,13 +1589,13 @@ 750 continue - do j=1,ntra - do i=minorig+1,nl - do il=1,ncum - if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then - traent(il,i,i,j)=tra(il,nk(il),j) - endif - enddo - enddo - enddo +c do j=1,ntra +c do i=minorig+1,nl +c do il=1,ncum +c if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then +c traent(il,i,i,j)=tra(il,nk(il),j) +c endif +c enddo +c enddo +c enddo do 100 j=minorig,nl @@ -1764,13 +1763,12 @@ enddo ! il - do j=1,ntra - do il=1,ncum - if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) - : .and. csum(il,i).lt.m(il,i) ) then - traent(il,i,i,j)=tra(il,nk(il),j) - endif - enddo - enddo - +c do j=1,ntra +c do il=1,ncum +c if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) +c : .and. csum(il,i).lt.m(il,i) ) then +c traent(il,i,i,j)=tra(il,nk(il),j) +c endif +c enddo +c enddo 789 continue c @@ -1869,11 +1867,11 @@ enddo - do k=1,ntra - do i=1,nd - do il=1,ncum - trap(il,i,k)=tra(il,i,k) - enddo - enddo - enddo +c do k=1,ntra +c do i=1,nd +c do il=1,ncum +c trap(il,i,k)=tra(il,i,k) +c enddo +c enddo +c enddo c @@ -2103,10 +2101,10 @@ vp(il,i)=vp(il,i)/mp(il,i) - do j=1,ntra - trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1) +c do j=1,ntra +c trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1) ctestmaf : +trap(il,i,j)*(mp(il,i)-mp(il,i+1)) - : +tra(il,i,j)*(mp(il,i)-mp(il,i+1)) - trap(il,i,j)=trap(il,i,j)/mp(il,i) - end do +c : +tra(il,i,j)*(mp(il,i)-mp(il,i+1)) +c trap(il,i,j)=trap(il,i,j)/mp(il,i) +c end do else @@ -2125,7 +2123,7 @@ vp(il,i)=vp(il,i+1) - do j=1,ntra - trap(il,i,j)=trap(il,i+1,j) - end do +c do j=1,ntra +c trap(il,i,j)=trap(il,i+1,j) +c end do endif @@ -2226,11 +2224,11 @@ enddo - do j=1,ntra - do i=1,nd - do il=1,ncum - ftra(il,i,j)=0.0 - enddo - enddo - enddo +c do j=1,ntra +c do i=1,nd +c do il=1,ncum +c ftra(il,i,j)=0.0 +c enddo +c enddo +c enddo do i=1,nl @@ -2330,17 +2328,17 @@ enddo ! il - do j=1,ntra - do il=1,ncum - if (cvflag_grav) then - ftra(il,1,j)=ftra(il,1,j)+0.01*grav*work(il) - : *(mp(il,2)*(trap(il,2,j)-tra(il,1,j)) - : +am(il)*(tra(il,2,j)-tra(il,1,j))) - else - ftra(il,1,j)=ftra(il,1,j)+0.1*work(il) - : *(mp(il,2)*(trap(il,2,j)-tra(il,1,j)) - : +am(il)*(tra(il,2,j)-tra(il,1,j))) - endif - enddo - enddo +c do j=1,ntra +c do il=1,ncum +c if (cvflag_grav) then +c ftra(il,1,j)=ftra(il,1,j)+0.01*grav*work(il) +c : *(mp(il,2)*(trap(il,2,j)-tra(il,1,j)) +c : +am(il)*(tra(il,2,j)-tra(il,1,j))) +c else +c ftra(il,1,j)=ftra(il,1,j)+0.1*work(il) +c : *(mp(il,2)*(trap(il,2,j)-tra(il,1,j)) +c : +am(il)*(tra(il,2,j)-tra(il,1,j))) +c endif +c enddo +c enddo do j=2,nl @@ -2366,21 +2364,21 @@ enddo - do k=1,ntra - do j=2,nl - do il=1,ncum - if (j.le.inb(il)) then - - if (cvflag_grav) then - ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1) - : *(traent(il,j,1,k)-tra(il,1,k)) - else - ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1) - : *(traent(il,j,1,k)-tra(il,1,k)) - endif - - endif - enddo - enddo - enddo +c do k=1,ntra +c do j=2,nl +c do il=1,ncum +c if (j.le.inb(il)) then + +c if (cvflag_grav) then +c ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1) +c : *(traent(il,j,1,k)-tra(il,1,k)) +c else +c ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1) +c : *(traent(il,j,1,k)-tra(il,1,k)) +c endif + +c endif +c enddo +c enddo +c enddo c @@ -2488,21 +2486,21 @@ 1350 continue - do k=1,ntra - do il=1,ncum - if (i.le.inb(il)) then - dpinv=1.0/(ph(il,i)-ph(il,i+1)) - cpinv=1.0/cpn(il,i) - if (cvflag_grav) then - ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv - : *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k)) - : -ad(il)*(tra(il,i,k)-tra(il,i-1,k))) - else - ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv - : *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k)) - : -ad(il)*(tra(il,i,k)-tra(il,i-1,k))) - endif - endif - enddo - enddo +c do k=1,ntra +c do il=1,ncum +c if (i.le.inb(il)) then +c dpinv=1.0/(ph(il,i)-ph(il,i+1)) +c cpinv=1.0/cpn(il,i) +c if (cvflag_grav) then +c ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv +c : *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k)) +c : -ad(il)*(tra(il,i,k)-tra(il,i-1,k))) +c else +c ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv +c : *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k)) +c : -ad(il)*(tra(il,i,k)-tra(il,i-1,k))) +c endif +c endif +c enddo +c enddo do 480 k=1,i-1 @@ -2538,21 +2536,21 @@ 480 continue - do j=1,ntra - do k=1,i-1 - do il=1,ncum - if (i.le.inb(il)) then - dpinv=1.0/(ph(il,i)-ph(il,i+1)) - cpinv=1.0/cpn(il,i) - if (cvflag_grav) then - ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i) - : *(traent(il,k,i,j)-tra(il,i,j)) - else - ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i) - : *(traent(il,k,i,j)-tra(il,i,j)) - endif - endif - enddo - enddo - enddo +c do j=1,ntra +c do k=1,i-1 +c do il=1,ncum +c if (i.le.inb(il)) then +c dpinv=1.0/(ph(il,i)-ph(il,i+1)) +c cpinv=1.0/cpn(il,i) +c if (cvflag_grav) then +c ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i) +c : *(traent(il,k,i,j)-tra(il,i,j)) +c else +c ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i) +c : *(traent(il,k,i,j)-tra(il,i,j)) +c endif +c endif +c enddo +c enddo +c enddo do 490 k=i,nl+1 @@ -2581,21 +2579,21 @@ 490 continue - do j=1,ntra - do k=i,nl+1 - do il=1,ncum - if (i.le.inb(il) .and. k.le.inb(il)) then - dpinv=1.0/(ph(il,i)-ph(il,i+1)) - cpinv=1.0/cpn(il,i) - if (cvflag_grav) then - ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i) - : *(traent(il,k,i,j)-tra(il,i,j)) - else - ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i) - : *(traent(il,k,i,j)-tra(il,i,j)) - endif - endif ! i and k - enddo - enddo - enddo +c do j=1,ntra +c do k=i,nl+1 +c do il=1,ncum +c if (i.le.inb(il) .and. k.le.inb(il)) then +c dpinv=1.0/(ph(il,i)-ph(il,i+1)) +c cpinv=1.0/cpn(il,i) +c if (cvflag_grav) then +c ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i) +c : *(traent(il,k,i,j)-tra(il,i,j)) +c else +c ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i) +c : *(traent(il,k,i,j)-tra(il,i,j)) +c endif +c endif ! i and k +c enddo +c enddo +c enddo do 1400 il=1,ncum @@ -2654,23 +2652,22 @@ enddo - do j=1,ntra - do il=1,ncum - if (i.le.inb(il)) then - dpinv=1.0/(ph(il,i)-ph(il,i+1)) - cpinv=1.0/cpn(il,i) - - if (cvflag_grav) then - ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv - : *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j)) - : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j))) - else - ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv - : *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j)) - : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j))) - endif - endif ! i - enddo - enddo - +c do j=1,ntra +c do il=1,ncum +c if (i.le.inb(il)) then +c dpinv=1.0/(ph(il,i)-ph(il,i+1)) +c cpinv=1.0/cpn(il,i) + +c if (cvflag_grav) then +c ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv +c : *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j)) +c : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j))) +c else +c ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv +c : *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j)) +c : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j))) +c endif +c endif ! i +c enddo +c enddo 500 continue @@ -2715,15 +2712,15 @@ 503 continue - do j=1,ntra - do il=1,ncum - ex=0.1*ment(il,inb(il),inb(il)) - : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j)) - : /(ph(il,inb(il))-ph(il,inb(il)+1)) - ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex - ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j) - : +ex*(ph(il,inb(il))-ph(il,inb(il)+1)) - : /(ph(il,inb(il)-1)-ph(il,inb(il))) - enddo - enddo +c do j=1,ntra +c do il=1,ncum +c ex=0.1*ment(il,inb(il),inb(il)) +c : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j)) +c : /(ph(il,inb(il))-ph(il,inb(il)+1)) +c ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex +c ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j) +c : +ex*(ph(il,inb(il))-ph(il,inb(il)+1)) +c : /(ph(il,inb(il)-1)-ph(il,inb(il))) +c enddo +c enddo c @@ -2981,4 +2978,30 @@ end + SUBROUTINE cv3_tracer(nloc,len,ncum,nd,na, + & ment,sij,da,phi) + implicit none +c inputs: + integer ncum, nd, na, nloc,len + real ment(nloc,na,na),sij(nloc,na,na) +c ouputs: + real da(nloc,na),phi(nloc,na,na) +c local variables: + integer i,j,k +c + da(:,:)=0. +c + do j=1,na + do k=1,na + do i=1,ncum + da(i,j)=da(i,j)+(1.-sij(i,k,j))*ment(i,k,j) + phi(i,j,k)=sij(i,k,j)*ment(i,k,j) +c print *,'da',j,k,da(i,j),sij(i,k,j),ment(i,k,j) + end do + end do + end do + + return + end + SUBROUTINE cv3_uncompress(nloc,len,ncum,nd,ntra,idcum @@ -3051,11 +3074,11 @@ - do 2100 j=1,ntra - do 2110 k=1,nd ! oct3 - do 2120 i=1,ncum - ftra1(idcum(i),k,j)=ftra(i,k,j) - 2120 continue - 2110 continue - 2100 continue +c do 2100 j=1,ntra +c do 2110 k=1,nd ! oct3 +c do 2120 i=1,ncum +c ftra1(idcum(i),k,j)=ftra(i,k,j) +c 2120 continue +c 2110 continue +c 2100 continue return Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/orografi.F =================================================================== --- /LMDZ.3.3/branches/rel-LF/libf/phylmd/orografi.F (revision 494) +++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/orografi.F (revision 495) @@ -315,5 +315,5 @@ call gwprofil * ( nlon , nlev - * , kgwd , kdx + * , kgwd , kdx , ktest * , ikcrith, icrit * , paphm1, zrho , zstab , zvph @@ -343,6 +343,10 @@ c c - do 523 jl=1,kgwd - ji=kdx(jl) +c do 523 jl=1,kgwd +c ji=kdx(jl) +c Modif vectorisation 02/04/2004 + do 523 ji=kidia,kfdia + if(ktest(ji).eq.1) then + zdelp=paphm1(ji,jk+1)-paphm1(ji,jk) ztemp=-rg*(ztau(ji,jk+1)-ztau(ji,jk))/(zvph(ji,ilevp1)*zdelp) @@ -401,4 +405,5 @@ pte(ji,jk)=0.0 + endif 523 continue @@ -1007,5 +1012,5 @@ SUBROUTINE GWPROFIL * ( NLON, NLEV - * , kgwd, kdx + * , kgwd, kdx , ktest * , KKCRITH, KCRIT * , PAPHM1, PRHO , PSTAB , PVPH , PRI , PTAU @@ -1075,5 +1080,6 @@ integer nlon,nlev INTEGER KKCRITH(NLON),KCRIT(NLON) - * ,kdx(nlon) + * ,kdx(nlon) , ktest(nlon) + C REAL PAPHM1(NLON,NLEV+1), PSTAB(NLON,NLEV+1), @@ -1109,8 +1115,12 @@ ilevh=KLEV/3 C - DO 400 ji=1,kgwd - jl=kdx(ji) +c DO 400 ji=1,kgwd +c jl=kdx(ji) +c Modif vectorisation 02/04/2004 + DO 400 jl=kidia,kfdia + if (ktest(jl).eq.1) then Zoro(JL)=Psig(JL)*Pdmod(JL)/4./max(pvar(jl),1.0) ZTAU(JL,KLEV+1)=PTAU(JL,KLEV+1) + endif 400 CONTINUE @@ -1123,6 +1133,9 @@ 410 CONTINUE C - DO 411 ji=1,kgwd - jl=kdx(ji) +c DO 411 ji=1,kgwd +c jl=kdx(ji) +c Modif vectorisation 02/04/2004 + do 411 jl=kidia,kfdia + if (ktest(jl).eq.1) then IF(JK.GT.KKCRITH(JL)) THEN PTAU(JL,JK)=ZTAU(JL,KLEV+1) @@ -1132,4 +1145,5 @@ PTAU(JL,JK)=GRAHILO*ZTAU(JL,KLEV+1) ENDIF + endif 411 CONTINUE C @@ -1143,6 +1157,9 @@ 420 CONTINUE C - DO 421 ji=1,kgwd - jl=kdx(ji) +c DO 421 ji=1,kgwd +c jl=kdx(ji) +c Modif vectorisation 02/04/2004 + do 421 jl=kidia,kfdia + if(ktest(jl).eq.1) then IF(JK.LT.KKCRITH(JL)) THEN ZNORM(JL)=gkdrag*PRHO(JL,JK)*SQRT(PSTAB(JL,JK))*PVPH(JL,JK) @@ -1150,4 +1167,5 @@ ZDZ2(JL,JK)=PTAU(JL,JK+1)/max(ZNORM(JL),gssec) ENDIF + endif 421 CONTINUE C @@ -1157,6 +1175,10 @@ C - DO 431 ji=1,kgwd - jl=kdx(ji) +c DO 431 ji=1,kgwd +c jl=Kdx(ji) +c Modif vectorisation 02/04/2004 + do 431 jl=kidia,kfdia + if(ktest(jl).eq.1) then + IF(JK.LT.KKCRITH(JL)) THEN IF((PTAU(JL,JK+1).LT.GTSEC).OR.(JK.LE.KCRIT(JL))) THEN @@ -1178,4 +1200,5 @@ ENDIF ENDIF + endif 431 CONTINUE @@ -1185,14 +1208,22 @@ C REORGANISATION OF THE STRESS PROFILE AT LOW LEVEL - DO 530 ji=1,kgwd - jl=kdx(ji) +c DO 530 ji=1,kgwd +c jl=kdx(ji) +c Modif vectorisation 02/04/2004 + do 530 jl=kidia,kfdia + if(ktest(jl).eq.1) then ZTAU(JL,KKCRITH(JL))=PTAU(JL,KKCRITH(JL)) ZTAU(JL,NSTRA)=PTAU(JL,NSTRA) + endif 530 CONTINUE DO 531 JK=1,KLEV - DO 532 ji=1,kgwd - jl=kdx(ji) +c DO 532 ji=1,kgwd +c jl=kdx(ji) +c Modif vectorisation 02/04/2004 + do 532 jl=kidia,kfdia + if(ktest(jl).eq.1) then + IF(JK.GT.KKCRITH(JL))THEN @@ -1206,10 +1237,15 @@ ENDIF + endif 532 CONTINUE C REORGANISATION IN THE STRATOSPHERE - DO 533 ji=1,kgwd - jl=kdx(ji) +c DO 533 ji=1,kgwd +c jl=kdx(ji) +c Modif vectorisation 02/04/2004 + do 533 jl=kidia,kfdia + if(ktest(jl).eq.1) then + IF(JK.LT.NSTRA)THEN @@ -1221,10 +1257,15 @@ ENDIF + endif 533 CONTINUE C REORGANISATION IN THE TROPOSPHERE - DO 534 ji=1,kgwd - jl=kdx(ji) +c DO 534 ji=1,kgwd +c jl=kdx(ji) +c Modif vectorisation 02/04/2004 + do 534 jl=kidia,kfdia + if(ktest(jl).eq.1) then + IF(JK.LT.KKCRITH(JL).AND.JK.GT.NSTRA)THEN @@ -1236,4 +1277,5 @@ ENDIF + endif 534 CONTINUE