c c $Header$ c subroutine read_reanalyse(timestep,u,v,t,masse,ps,mode) c mode=0 variables naturelles c mode=1 variabels GCM IMPLICIT NONE #include "netcdf.inc" #include "dimensions.h" #include "paramet.h" #include "comgeom.h" #include "comvert.h" integer nlevnc parameter (nlevnc=15) C pour annee 2000 parameter (nlevnc=21) integer timestep,mode,l real u(iip1,jjp1,llm),v(iip1,jjm,llm) real t(iip1,jjp1,llm),ps(iip1,jjp1) real masse(iip1,jjp1,llm),pk(iip1,jjp1,llm) integer ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps integer varidpl save ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps save varidpl real*4 unc(iip1,jjp1,nlevnc),vnc(iip1,jjm,nlevnc) real*4 tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1) real*4 pl(nlevnc) integer start(4),count(4),status real rcode logical first,ncep save first,ncep cMod 11 2 99 data first,ncep/.true.,.false./ c data first,ncep/.true.,.true./ print*,'ok0' c ----------------------------------------------------------------- c Initialisation de la lecture des fichiers c ----------------------------------------------------------------- if (first) then ncidu=NCOPN('u.nc',NCNOWRIT,rcode) varidu=NCVID(ncidu,'UWND',rcode) print*,'ncidu,varidu',ncidu,varidu if (ncep) then varidpl=NCVID(ncidu,'LEVEL',rcode) else varidpl=NCVID(ncidu,'PRESSURE',rcode) endif print*,'ncidu,varidpl',ncidu,varidpl ncidv=NCOPN('v.nc',NCNOWRIT,rcode) varidv=NCVID(ncidv,'VWND',rcode) print*,'ncidv,varidv',ncidv,varidv ncidt=NCOPN('T.nc',NCNOWRIT,rcode) varidt=NCVID(ncidt,'AIR',rcode) print*,'ncidt,varidt',ncidt,varidt c ncidps=NCOPN('ps.nc',NCNOWRIT,rcode) c varidps=NCVID(ncidps,'SP',rcode) c print*,'ncidps,varidps',ncidps,varidps endif print*,'ok1' c ----------------------------------------------------------------- c lecture des champs u, v, T, ps c ----------------------------------------------------------------- c niveaux de pression print*,'WARNING!!! Il n y a pas de test de coherence' print*,'sur le nombre de niveaux verticaux dans le fichier nc' #ifdef NC_DOUBLE status=NF_GET_VARA_DOUBLE(ncidu,varidpl,1,nlevnc,pl) #else status=NF_GET_VARA_REAL(ncidu,varidpl,1,nlevnc,pl) #endif start(1)=1 start(2)=1 start(3)=1 start(4)=timestep count(1)=iip1 count(2)=jjp1 count(3)=nlevnc count(4)=1 c u #ifdef NC_DOUBLE status=NF_GET_VARA_DOUBLE(ncidu,varidu,start,count,unc) #else status=NF_GET_VARA_REAL(ncidu,varidu,start,count,unc) #endif print*,'WARNING!!! Correction bidon pour palier a un ' print*,'probleme dans la creation des fichiers nc' call correctbid(iim,jjp1*nlevnc,unc) c call dump2d(iip1,jjp1,unc,'UNC COUCHE 1 ') c T #ifdef NC_DOUBLE status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,tnc) #else status=NF_GET_VARA_REAL(ncidt,varidt,start,count,tnc) #endif call correctbid(iim,jjp1*nlevnc,tnc) c call dump2d(iip1,jjp1,tnc,'TNC COUCHE 1 ') count(2)=jjm c v #ifdef NC_DOUBLE status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,vnc) #else status=NF_GET_VARA_REAL(ncidv,varidv,start,count,vnc) #endif call correctbid(iim,jjm*nlevnc,vnc) c call dump2d(iip1,jjm,vnc,'VNC COUCHE 1 ') start(3)=timestep start(4)=0 count(2)=jjp1 count(3)=1 count(4)=0 c ps C#ifdef NC_DOUBLE c status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,psnc) C#else c status=NF_GET_VARA_REAL(ncidps,varidps,start,count,psnc) C#endif c call dump2d(iip1,jjp1,psnc,'PSNC COUCHE 1 ') c call correctbid(iim,jjp1,psnc) c Transformations do l=1,nlevnc pl(l)=pl(l)*100. enddo if(first) then do l=1,nlevnc print*,'PL(',l,')=',pl(l) enddo endif call reanalyse2nat(nlevnc,unc,vnc,tnc,psnc,pl,u,v,t s ,ps,masse,pk) c call dump2d(iip1,jjm,v,'V COUCHE 1 ') if(mode.eq.1) call nat2gcm(u,v,t,pk,u,v,t) print*,'TIMESTEP ',timestep if(mode.ne.1) stop'mode pas egal 0' c call dump2d(iip1,jjm,v,'VCOV COUCHE 1 ') c Lignes introduites a une epoque pour un probleme oublie... c do l=1,llm c do i=1,iip1 c v(i,1,l)=0. c v(i,jjm,l)=0. c enddo c enddo first=.false. return end c ----------------------------------------------------------------- subroutine reanalyse2nat(nlevnc,unc,vnc,tnc,psnc,pl,u,v,t s ,ps,masse,pk) c ----------------------------------------------------------------- c Inversion Nord/sud de la grille + interpollation sur les niveaux c verticaux du modele. c ----------------------------------------------------------------- implicit none #include "dimensions.h" #include "paramet.h" #include "comgeom2.h" #include "comvert.h" #include "comconst.h" integer nlevnc real u(iip1,jjp1,llm),v(iip1,jjm,llm) real t(iip1,jjp1,llm),ps(iip1,jjp1) real pl(nlevnc) real unc(iip1,jjp1,nlevnc),vnc(iip1,jjm,nlevnc) real tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1) real zu(iip1,jjp1,llm),zv(iip1,jjm,llm) real zt(iip1,jjp1,llm) real pext(iip1,jjp1,llm) real pbarx(iip1,jjp1,llm),pbary(iip1,jjm,llm) real plunc(iip1,jjp1,llm),plvnc(iip1,jjm,llm) real plsnc(iip1,jjp1,llm) REAL alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm) real p(iip1,jjp1,llmp1),pk(iip1,jjp1,llm),pks(iip1,jjp1) real pkf(iip1,jjp1,llm) real masse(iip1,jjp1,llm),pls(iip1,jjp1,llm) real prefkap,unskap integer i,j,l c ----------------------------------------------------------------- c calcul de la pression au milieu des couches. c ----------------------------------------------------------------- c do j=1,jjp1 c do i=1,iim c ps(i,j)=psnc(i,jjp1+1-j) c enddo c ps(iip1,j)=ps(1,j) c enddo CALL pression( ip1jmp1, ap, bp, ps, p ) call massdair(p,masse) CALL exner_hyb(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf) c .... Calcul de pls , pression au milieu des couches ,en Pascals unskap=1./kappa prefkap = preff ** kappa c PRINT *,' Pref kappa unskap ',preff,kappa,unskap DO l = 1, llm DO j=1,jjp1 DO i =1, iip1 pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap ENDDO ENDDO ENDDO c ----------------------------------------------------------------- c calcul des pressions pour les grilles u et v c ----------------------------------------------------------------- do l=1,llm do j=1,jjp1 do i=1,iip1 pext(i,j,l)=pls(i,j,l)*aire(i,j) enddo enddo enddo call massbar(pext, pbarx, pbary ) do l=1,llm do j=1,jjp1 do i=1,iip1 plunc(i,jjp1+1-j,l)=pbarx(i,j,l)/aireu(i,j) plsnc(i,jjp1+1-j,l)=pls(i,j,l) enddo enddo enddo do l=1,llm do j=1,jjm do i=1,iip1 plvnc(i,jjm+1-j,l)=pbary(i,j,l)/airev(i,j) enddo enddo enddo c ----------------------------------------------------------------- call pres2lev(unc,zu,nlevnc,llm,pl,plunc,iip1,jjp1) call pres2lev(vnc,zv,nlevnc,llm,pl,plvnc,iip1,jjm ) call pres2lev(tnc,zt,nlevnc,llm,pl,plsnc,iip1,jjp1) c call dump2d(iip1,jjp1,ps,'PS ') c call dump2d(iip1,jjp1,psu,'PS ') c call dump2d(iip1,jjm,psv,'PS ') c Inversion Nord/Sud do l=1,llm do j=1,jjp1 do i=1,iim u(i,j,l)=zu(i,jjp1+1-j,l) t(i,j,l)=zt(i,jjp1+1-j,l) enddo u(iip1,j,l)=u(1,j,l) t(iip1,j,l)=t(1,j,l) enddo enddo do l=1,llm do j=1,jjm do i=1,iim v(i,j,l)=zv(i,jjm+1-j,l) enddo v(iip1,j,l)=v(1,j,l) enddo enddo return end subroutine nat2gcm(u,v,t,pk,ucov,vcov,teta) implicit none #include "dimensions.h" #include "paramet.h" #include "comgeom2.h" #include "comconst.h" #include "comvert.h" real u(iip1,jjp1,llm),v(iip1,jjm,llm) real t(iip1,jjp1,llm),pk(iip1,jjp1,llm) real ucov(iip1,jjp1,llm),vcov(iip1,jjm,llm) real teta(iip1,jjp1,llm) integer i,j,l c calcul de ucov et de la temperature potentielle do l=1,llm do j=1,jjp1 do i=1,iim ucov(i,j,l)=u(i,j,l)*cu(i,j) teta(i,j,l)=t(i,j,l)*cpp/pk(i,j,l) enddo ucov(iip1,j,l)=ucov(1,j,l) teta(iip1,j,l)=teta(1,j,l) enddo do i=1,iip1 ucov(i,1,l)=0. ucov(i,jjp1,l)=0. teta(i,1,l)=teta(1,1,l) teta(i,jjp1,l)=teta(1,jjp1,l) enddo enddo c calcul de ucov do l=1,llm do j=1,jjm do i=1,iim vcov(i,j,l)=v(i,j,l)*cv(i,j) enddo vcov(iip1,j,l)=vcov(1,j,l) enddo enddo c call dump2d(iip1,jjp1,teta,'TETA EN BAS ') c call dump2d(iip1,jjp1,teta(1,1,llm),'TETA EN HAUT ') return end subroutine correctbid(iim,nl,x) integer iim,nl real x(iim+1,nl) integer i,l real zz do l=1,nl do i=2,iim-1 if(abs(x(i,l)).gt.1.e10) then zz=0.5*(x(i-1,l)+x(i+1,l)) c print*,'correction ',i,l,x(i,l),zz x(i,l)=zz endif enddo enddo return end