SUBROUTINE redecoupe(irec,massemn,pbarun,pbarvn,wn,tetan,phin, s nrec,avant,airefi, s zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz, s yu1,yv1,ftsol,pctsrf, s frac_impa,frac_nucl,phisn) IMPLICIT NONE #include "dimensions.h" #include "paramet.h" #include "comvert.h" #include "comconst.h" #include "comgeom2.h" #include "tracstoke.h" integer irec,nrec,i,j integer ig,l integer imo,jmo,imn,jmn,ii,jj,ig parameter (imn=iim,jmn=jjm,imo=imn/2,jmo=(jmn+1)/2) integer ngrido,ngridn parameter(ngrido=(jmo-1)*imo+2,ngridn=(jmn-1)*imn+2) INTEGER nbsrf PARAMETER (nbsrf=4) ! nombre de sous-fractions pour une maille real zmfd(ngridn,llm),zde_d(ngridn,llm),zen_d(ngridn,llm) real zmfu(ngridn,llm),zde_u(ngridn,llm),zen_u(ngridn,llm) real coefkz(ngridn,llm) real frac_impa(ngridn,llm),frac_nucl(ngridn,llm) real yu1(ngridn), yv1(ngridn) real ftsol(ngridn,nbsrf),pctsrf(ngridn,nbsrf) integer imfu,imfd,ien_u,ide_u, s ien_d,ide_d, s icoefkz,izu1,izv1, s itsol,ipsf, s ilei, ilec parameter(imfu=1,imfd=llm+1,ien_u=2*llm+1,ide_u=3*llm+1, s ien_d=4*llm+1,ide_d=5*llm+1, s icoefkz=6*llm+1, s ilei=7*llm+1,ilec=8*llm+1, s izu1=9*llm+1,izv1=9*llm+2, s itsol=9*llm+3,ipsf=9*llm+3+nbsrf) logical avant real massefi(ngridn,llm) real massemn(imn+1,jmn+1,llm),tetan(imn+1,jmn+1,llm) real pbarun(imn+1,jmn+1,llm),pbarvn(imn+1,jmn,llm) real wn(imn+1,jmn+1,llm),phin(imn+1,jmn+1,llm) real phisn(imn+1,jmn+1) real massemo(imo+1,jmo+1,llm),tetao(imo+1,jmo+1,llm) real pbaruo(imo+1,jmo+1,llm),pbarvo(imo+1,jmo,llm) real wo(imo+1,jmo+1,llm),phio(imo+1,jmo+1,llm) real phiso(imo+1,jmo+1) real pbarvst(imo+1,jmo+1,llm) real airefi(ngridn) real xlecn(ngridn,9*llm+2+2*nbsrf),tmpn(imn+1,jmn+1) real xleco(ngrido,9*llm+2+2*nbsrf),tmpo(imo+1,jmo+1) real zcontrole(ngridn),zmass,tmpdyn(imn+1,jmn+1),zflux real ziadvtrac,zrec,ziadvtrac2,zrec2 real zim,zjm,zlm,zklon,zklev real zpi c longitudes et latitudes lues real rlonul(1:imo+1),rlatvl(1:jmo) real rlonvl(1:imo+1),rlatul(1:jmo+1) c longitudes et latitudes anciennes real rlonuo(0:imo+1),rlatvo(0:jmo+1) c longitudes et latitudes nouvelles real rlonun(0:imn+1),rlatvn(0:jmn+1) real aireo(imo+1,jmo+1) integer ndecx(imo+1),ndecy(jmo+1) real alphax(imn+1),alphay(jmn+1) real alphaxo(imo+1) real alpha(imn+1,jmn+1) real aa,uu(0:imo+1),vv(imo+1,0:jmo+1) integer iest(imo+1),iouest(imo+1) integer jsud(jmo+1),jnord(jmo+1) integer in,io,jn,jo,l,iin,jjn integer i,j real dlatm,dlatp,dlonm,dlonp zpi=2.*asin(1.) c================================================================== c Si le numero du record est 0 alors: INITIALISATION c================================================================== c print*,'ENTREE DANS LECTFLUX' print*,'IREC=',IREC if(irec.eq.0) then print*,'IREC==',0 C test call inigeom c================================================================== c Definition des surdecoupages dans les deux directions c================================================================== ndecx(1)=1 do io=2,imo ndecx(io)=2 enddo ndecx(imo+1)=1 ndecy(1)=1 do jo=2,jmo ndecy(jo)=2 enddo ndecy(jmo+1)=1 ii=0 do io=1,imo+1 ii=ii+ndecx(io) enddo if(ii.ne.iim) then print*,'ii=',ii,' iim=',iim stop endif jj=0 do jo=1,jmo+1 jj=jj+ndecy(jo) enddo if(jj.ne.jjp1) then print*,'jj=',jj,' jjm=',jjm stop endif c================================================================== c Calcul des jsud,... correspondant aux intersections des c grilles. c================================================================== iest(1)=0 do io=2,imo+1 iest(io)=iest(io-1)+ndecx(io-1) iouest(io-1)=iest(io) enddo iouest(imo+1)=iest(imo+1)+ndecx(imo+1) jnord(1)=0 do jo=2,jmo+1 jnord(jo)=jnord(jo-1)+ndecy(jo-1) jsud(jo-1)=jnord(jo) enddo jsud(jmo+1)=jnord(jmo+1)+ndecy(jmo+1) c================================================================== c ouverture des fichiers, lecture des entetes c================================================================== c Fichier fluxmass #ifdef CRAY CALL ASSIGN("assign -N ieee -F null f:fluxmass") #endif open(47,file='fluxmass',form='unformatted', s access='direct' s ,recl=4*6*(imo+1)*(jmo+1)*llm) read(47,rec=1) zrec,dtvr,ziadvtrac,zim,zjm,zlm s ,rlonul,rlonvl,rlatul,rlatvl,aireo s ,phiso print*,'zrec,dtvr,ziadvtrac,zim,zjm,zlm' print*,zrec,dtvr,ziadvtrac,zim,zjm,zlm if((imo-nint(zim))*(jmo-nint(zjm)).ne.0) then print*,'Modifier les dimensions dans redecoupe ' print*,'Mettre imo=',zim,' jmo=',zjm stop endif c Fichier physique c Fichier lessivage (supprime les donnees utiles sont dans "physique") #ifdef CRAY CALL ASSIGN("assign -N ieee -F null f:physique") #endif open(49,file='physique',form='unformatted', s access='direct' s ,recl=4*ngrido*(9*llm+2+2*nbsrf)) read(49,rec=1) zrec2,ziadvtrac2,zklon,zklev print*,'Entete du fichier physique' print*,zrec2,ziadvtrac2,zklon,zklev nrec=zrec istdyn=ziadvtrac istphy=ziadvtrac2 c================================================================== c Definition des anciennes latitudes et longitudes c (qui pourraient etre relues plus tard) c================================================================== rlonuo(0)=-zpi do io=1,imo c rlonuo(io)=2.*zpi/FLOAT(imo)*(io+0.5-0.5*FLOAT(imo)-1.) c print*,'LON ',io,rlonuo(io),rlonul(io) rlonuo(io)=rlonul(io) enddo rlonuo(imo+1)=zpi rlatvo(0)=zpi/2. do jo=1,jmo c rlatvo(jo)=zpi/FLOAT(jmo)*(0.5*FLOAT(jmo)+1.-jo-0.5) c print*,'LAT ',jo,rlatvo(jo),rlatvl(jo) rlatvo(jo)=rlatvl(jo) enddo rlatvo(jmo+1)=-zpi/2. c do jo=1,jmo+1 c do io=1,imo+1 c aireo(io,jo)=rad*rad c s *(rlonuo(io)-rlonuo(io-1)) c s *(sin(rlatvo(jo-1))-sin(rlatvo(jo))) c aireo(io,jo)=airel(io,jo) c enddo c aireo(1,jo)=aireo(1,jo)+aireo(imo+1,jo) c aireo(imo+1,jo)=aireo(1,jo) c enddo do io=2,imo alphaxo(io)=1. enddo alphaxo(1)=(rlonuo(1)-rlonuo(0)) s /(rlonuo(1)-rlonuo(0)+rlonuo(imo+1)-rlonuo(imo)) alphaxo(imo+1)=1.-alphaxo(1) c================================================================== c Definition des nouvelles latitudes et longitudes c================================================================== rlonun(0)=-zpi do io=1,imo+1 do iin=1,iouest(io)-iest(io) in=iin+iest(io) rlonun(in)= s rlonuo(io-1)+iin*(rlonuo(io)-rlonuo(io-1)) s /ndecx(io) alphax(in)=alphaxo(io)/ndecx(io) print787,io,rlonuo(io-1)*180./zpi,in s ,iest(io),iouest(io),rlonun(in)*180./zpi,alphax(in) enddo enddo rlatvn(0)=0.5*zpi do jo=1,jmo+1 print*,'jo=',jo do jjn=1,jsud(jo)-jnord(jo) jn=jnord(jo)+jjn rlatvn(jn)=rlatvo(jo-1)+jjn*(rlatvo(jo)-rlatvo(jo-1)) s /ndecy(jo) alphay(jn)=(sin(rlatvn(jn-1))-sin(rlatvn(jn))) s /(sin(rlatvo(jo-1))-sin(rlatvo(jo))) print787,jo,rlatvo(jo-1)*180./zpi,jn s ,jnord(jo),jsud(jo),rlatvn(jn)*180./zpi,alphay(jn) enddo enddo 787 format(i5,f10.2,3(i5),2(f10.2)) do in=1,imn rlonu(in)=rlonun(in) rlonv(in)=0.5*(rlonun(in)+rlonun(in-1)) enddo rlonv(imn+1)=rlonv(1)+2.*zpi rlonu(imn+1)=rlonu(1)+2.*zpi do jn=1,jmn rlatv(jn)=rlatvn(jn) enddo do jn=1,jmn+1 rlatu(jn)=0.5*(rlatvn(jn-1)+rlatvn(jn)) enddo do jn=1,jmn+1 do in=1,imn alpha(in,jn)=alphax(in)*alphay(jn) enddo alpha(imn+1,jn)=0. enddo c call dump2d(iip1,jjp1,alpha,'ALPHA ') c . on a : cu(i,j) = rad * COS(y) * dx/dX . c . cv( j ) = rad * dy/dY . c A 1 point scalaire P (i,j) de la grille, reguliere en (X,Y) , sont c affectees 4 aires entourant P , calculees respectivement aux points c ( i + 1/4, j - 1/4 ) : aireij1 (i,j) c ( i + 1/4, j + 1/4 ) : aireij2 (i,j) c ( i - 1/4, j + 1/4 ) : aireij3 (i,j) c ( i - 1/4, j - 1/4 ) : aireij4 (i,j) c c . V c c aireij4 . . aireij1 c c U . . P . U c c aireij3 . . aireij2 c c . V do j=1,jjp1 do i=1,iim dlonp=rlonun(i)-rlonv(i) dlonm=rlonv(i)-rlonun(i-1) dlatp=sin(rlatvn(j-1))-sin(rlatu(j)) dlatm=sin(rlatu(j))-sin(rlatvn(j)) aireij1 ( i,j ) = rad*rad*dlatp*dlonp aireij2 ( i,j ) = rad*rad*dlatm*dlonp aireij3 ( i,j ) = rad*rad*dlatm*dlonm aireij4 ( i,j ) = rad*rad*dlatp*dlonm aire ( i,j ) = aireij1(i,j) + aireij2(i,j) + aireij3(i,j) + * aireij4(i,j) alpha1 ( i,j ) = aireij1(i,j) / aire(i,j) alpha2 ( i,j ) = aireij2(i,j) / aire(i,j) alpha3 ( i,j ) = aireij3(i,j) / aire(i,j) alpha4 ( i,j ) = aireij4(i,j) / aire(i,j) alpha1p2( i,j ) = alpha1 (i,j) + alpha2 (i,j) alpha1p4( i,j ) = alpha1 (i,j) + alpha4 (i,j) alpha2p3( i,j ) = alpha2 (i,j) + alpha3 (i,j) alpha3p4( i,j ) = alpha3 (i,j) + alpha4 (i,j) enddo aireij1(iip1,j)=aireij1(1,j) aireij2(iip1,j)=aireij2(1,j) aireij3(iip1,j)=aireij3(1,j) aireij4(iip1,j)=aireij4(1,j) aire(iip1,j)=aire(1,j) alpha1(iip1,j)=alpha1(1,j) alpha2(iip1,j)=alpha2(1,j) alpha3(iip1,j)=alpha3(1,j) alpha4(iip1,j)=alpha4(1,j) alpha1p2(iip1,j)=alpha1p2(1,j) alpha1p4(iip1,j)=alpha1p4(1,j) alpha2p3(iip1,j)=alpha2p3(1,j) alpha3p4(iip1,j)=alpha3p4(1,j) enddo c call dump2d(iip1,jjp1,aire,'AIRE ') c do jn=1,jjp1 c do in=1,iim c aire(in,jn)=rad*rad*(sin(rlatvn(jn-1))-sin(rlatvn(jn))) c s *(rlonun(in)-rlonun(in-1)) c unsaire(in,jn)=1./aire(in,jn) c enddo c aire(iip1,jn)=aire(1,jn) c unsaire(iip1,jn)=unsaire(1,jn) c enddo c call dump2d(iip1,jjp1,aire,'AIRE2 ') DO 42 j = 1,jjp1 DO 41 i = 1,iim unsaire(i,j) = 1./ aire(i,j) aireu (i,j) = aireij1(i,j) + aireij2(i,j) + aireij4(i+1,j) + * aireij3(i+1,j) 41 CONTINUE aireu (iip1,j) = aireu (1,j) unsaire(iip1,j) = unsaire(1,j) 42 CONTINUE DO 48 j = 1,jjm DO i=1,iim airev(i,j) = aireij2(i,j)+ aireij3(i,j)+ aireij1(i,j+1) + * aireij4(i,j+1) ENDDO airev (iip1,j) = airev(1,j) 48 CONTINUE apoln=0. apols=0. do i=1,iim apoln=apoln+aire(i,1) apols=apols+aire(i,jjp1) enddo do jn=1,jjp1 do in=1,iim cu(in,jn)=rad*cos(rlatu(jn))*(rlonv(in+1)-rlonv(in)) enddo cu(iip1,jn)=cu(1,jn) enddo do jn=1,jjm do in=1,iim+1 cv(in,jn)=rad*(rlatu(jn)-rlatu(jn+1)) enddo enddo c================================================================== c Fin des initialisations else ! irec=0 c================================================================== c----------------------------------------------------------------------- c Lecture des fichiers fluxmass et physique: c ----------------------------------------------------- c================================================================== c Variables dynamiques c================================================================== read(47,rec=irec) massemo,pbaruo,pbarvst,wo,tetao,phio do l=1,llm do j=1,jmo do i=1,imo+1 pbarvo(i,j,l)=pbarvst(i,j,l) enddo enddo enddo do l=1,llm do jo=1,jmo+1 do io=1,imo+1 do jn=jnord(jo)+1,jsud(jo) do in=iest(io)+1,iouest(io) wn(in,jn,l)=alpha(in,jn)*wo(io,jo,l) massemn(in,jn,l)=alpha(in,jn) s *massemo(io,jo,l) tetan(in,jn,l)=tetao(io,jo,l) phin(in,jn,l)=phio(io,jo,l) enddo enddo enddo enddo do jn=1,jmn+1 wn(imn+1,jn,l)=wn(1,jn,l) massemn(imn+1,jn,l)=massemn(1,jn,l) tetan(imn+1,jn,l)=tetan(1,jn,l) phin(imn+1,jn,l)=phin(1,jn,l) enddo enddo do l=1,llm do jo=1,jmo+1 uu(imo+1)=0.5*(pbaruo(imo,jo,l)+pbaruo(imo+1,jo,l)) uu(0)=uu(imo+1) do io=1,imo uu(io)=pbaruo(io,jo,l) enddo do io=1,imo+1 do jn=jnord(jo)+1,jsud(jo) aa=0. do in=iest(io)+1,iouest(io) aa=aa+alphax(in) pbarun(in,jn,l)=alphay(jn)* s (uu(io-1)+aa*(uu(io)-uu(io-1))) enddo enddo enddo enddo do jn=1,jmn+1 pbarun(imn+1,jn,l)=pbarun(1,jn,l) enddo enddo do l=1,llm do jo=1,jmo do io=1,imo+1 vv(io,jo)=pbarvo(io,jo,l) enddo enddo do io=1,imo+1 vv(io,0)=0. vv(io,jmo+1)=0. enddo do jo=1,jmo+1 do io=1,imo+1 aa=0. c do jn=jnord(jo)+1,max(jsud(jo),jmo) do jn=jnord(jo)+1,min(jsud(jo),jmn) aa=aa+alphay(jn) do in=iest(io)+1,iouest(io) pbarvn(in,jn,l)=alphax(in)* s (vv(io,jo-1)+aa*(vv(io,jo)-vv(io,jo-1))) enddo enddo enddo enddo do jn=1,jmn pbarvn(iip1,jn,l)=pbarvn(1,jn,l) enddo enddo c================================================================== c Variables physiques c================================================================== read(49,rec=irec) ((xleco(ig,l),ig=1,ngrido), s l=1,9*llm+2+2*nbsrf) c================================================================== c Passage a la nouvelle grille c================================================================== do l=1,9*llm+2+2*nbsrf c passage aa la grille dynamique ancienne do io=1,imo+1 tmpo(io,1)=xleco(1,l) tmpo(io,jmo+1)=xleco(ngrido,l) enddo do jo=2,jmo do io=1,imo tmpo(io,jo)=xleco((jo-2)*imo+io+1,l) enddo tmpo(imo+1,jo)=tmpo(1,jo) enddo c passage a la grillle dynamique nouvelle do jo=1,jmo+1 do io=1,imo+1 do jn=jnord(jo)+1,jsud(jo) do in=iest(io)+1,iouest(io) tmpn(in,jn)=tmpo(io,jo) enddo enddo enddo enddo c passage a la grille physique nouvelle xlecn(1,l)=tmpn(1,1) xlecn(ngridn,l)=tmpn(1,jmn+1) do jn=2,jmn do in=1,imn xlecn((jn-2)*imn+in+1,l)=tmpn(in,jn) enddo enddo enddo c================================================================== do l=1,llm do ig=1,ngridn coefkz(ig,l)=xlecn(ig,icoefkz+l-1) frac_impa(ig,l)=xlecn(ig,ilei+l-1) frac_nucl(ig,l)=xlecn(ig,ilec+l-1) enddo enddo do l=1,nbsrf do ig=1,ngridn ftsol(ig,l)=xlecn(ig,itsol+l-1) pctsrf(ig,l)=xlecn(ig,ipsf+l-1) enddo enddo do ig=1,ngridn yv1(ig)=xlecn(ig,izv1) yu1(ig)=xlecn(ig,izu1) enddo C if(avant) then c Simu directe do l=1,llm do ig=1,ngridn zmfu(ig,l)=xlecn(ig,imfu+l-1) zmfd(ig,l)=xlecn(ig,imfd+l-1) zde_u(ig,l)=xlecn(ig,ide_u+l-1) zen_u(ig,l)=xlecn(ig,ien_u+l-1) zde_d(ig,l)=xlecn(ig,ide_d+l-1) zen_d(ig,l)=xlecn(ig,ien_d+l-1) enddo enddo else c Simu retro do l=1,llm do ig=1,ngridn zmfd(ig,l)=-xlecn(ig,imfu+l-1) zmfu(ig,l)=-xlecn(ig,imfd+l-1) zen_d(ig,l)=xlecn(ig,ide_u+l-1) zde_d(ig,l)=xlecn(ig,ien_u+l-1) zen_u(ig,l)=xlecn(ig,ide_d+l-1) zde_u(ig,l)=xlecn(ig,ien_d+l-1) enddo enddo endif c----------------------------------------------------------------------- c PETIT CONTROLE SUR LES FLUX CONVECTIFS... c----------------------------------------------------------------------- call gr_dyn_fi(llm,iip1,jjp1,ngridn,massemn,massefi) print*,'Ap redec irec' do ig=1,ngridn zcontrole(ig)=1. enddo c zmass=(max(massemn(ig,l),massemn(ig,l-1))/airefi(ig) do l=2,llm do ig=1,ngridn zmass=max(massefi(ig,l),massefi(ig,l-1))/airefi(ig) zflux=max(abs(zmfu(ig,l)),abs(zmfd(ig,l)))*dtphys if(zflux.gt.0.9*zmass) then zcontrole(ig)=min(zcontrole(ig),0.9*zmass/zflux) endif enddo enddo do ig=1,ngridn if(zcontrole(ig).lt.0.99999) then print*,'ATTENTION !!! on reduit les flux de masse ' print*,'convectifs au point ig=',ig endif enddo call gr_fi_dyn(1,ngridn,iip1,jjp1,zcontrole,tmpdyn) do l=1,llm do ig=1,ngridn zmfu(ig,l)=zmfu(ig,l)*zcontrole(ig) zmfd(ig,l)=zmfd(ig,l)*zcontrole(ig) zen_u(ig,l)=zen_u(ig,l)*zcontrole(ig) zde_u(ig,l)=zde_u(ig,l)*zcontrole(ig) zen_d(ig,l)=zen_d(ig,l)*zcontrole(ig) zde_d(ig,l)=zde_d(ig,l)*zcontrole(ig) enddo enddo endif ! irec=0 RETURN END