c c $Header$ c PROGRAM offlinenc USE ioipsl IMPLICIT NONE #include "dimensions.h" #include "paramet.h" #include "comconst.h" #include "comdissnew.h" #include "comvert.h" #include "comgeom2.h" #include "logic.h" #include "temps.h" #include "control.h" #include "ener.h" #include "description.h" #include "serre.h" #include "tracstoke.h" integer ngridmx,nbsrf PARAMETER (ngridmx=iim*(jjm-1)+2,nbsrf=4) c----------------------------------------------------------------------- INTEGER*4 iday ! jour julien REAL time ! Heure de la journee en fraction d'1 jour REAL dtdyn,dtav INTEGER iday_step INTEGER longcles PARAMETER ( longcles = 20 ) REAL clesphy0( longcles ) logical clefinit c variables dynamiques REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants REAL q(iip1,jjp1,llm,nqmx) ! champs advectes real qmoy(iip1,jjp1,nqmx) REAL pks(iip1,jjp1) ! exner (f pour filtre) REAL phis(iip1,jjp1) ! geopotentiel au sol REAL phi(iip1,jjp1,llm) ! geopotentiel REAL w(iip1,jjp1,llm) ! vitesse verticale C C variables physiques real zmfd(ngridmx,llm),zde_d(ngridmx,llm),zen_d(ngridmx,llm) REAL qfi(ngridmx,llm,nqmx) ! champs advectes real zmfu(ngridmx,llm),zde_u(ngridmx,llm),zen_u(ngridmx,llm) real t(ngridmx,llm) real coefkz(ngridmx,llm) real frac_impa(ngridmx,llm),frac_nucl(ngridmx,llm) real pphis (ngridmx) REAL airefi(ngridmx) REAL latfi(ngridmx),lonfi(ngridmx) ! latitude et long pour chaque point real yu1(ngridmx), yv1(ngridmx) real ftsol(ngridmx,nbsrf),pctsrf(ngridmx,nbsrf) C character*10 file c variables dynamiques intermediaire pour le transport REAL pbaru(iip1,jjp1,llm),pbarv(iip1,jjm,llm) !flux de masse REAL masse(iip1,jjp1,llm) REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm) REAL teta(iip1,jjp1,llm) INTEGER itau,irec,nrec,ldec,itauav integer recmin,recmax,recstkmin,recstkmax character*1 injecte,stoke character*10 nom real tinj2,tinj1 integer iiinj integer jour0,isplit,nsplit_dyn,nsplit_phy,nsplit logical debut,lectstart,rnpb,lafin EXTERNAL inidissip,iniconst,inifilr EXTERNAL inigeom EXTERNAL exner_hyb,pression EXTERNAL defrun_new REAL time_0 ! temps de debut de simulation (% journee) character*2 str2 INTEGER iq,j,i,l,nq,mode real paprs(ngridmx,llm+1),pplay(ngridmx,llm),tempfi(ngridmx,llm) real p(iip1,jjp1,llmp1),pk(iip1,jjp1,llm) real lon_stat(nqmx),lat_stat(nqmx),zlon,zlat integer i_stat(nqmx),j_stat(nqmx) real intensite(nqmx) integer iqmin,iqmax integer npasinject data npasinject/1/ integer histid, histvid, histaveid real time_step, t_wrt, t_ops integer itime_step character*80 dynhist_file, dynhistave_file logical avant logical debutphy,redecoupage DATA debutphy/.true./ real unskap,pksurcp,smass(iip1,jjp1) integer ig0,ii,jj,iii,ij REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm) REAL ps(iip1,jjp1),pkf(ip1jmp1,llm) integer ntau, iinjec2,iij real zttt c----------------------------------------------------------------------- c Initialisations: c ---------------- descript = 'Run GCM LMDZ' tinj1=0. tinj2=0. rad = 6400000 omeg = 7.272205e-05 g = 9.8 kappa = 0.285716 daysec = 86400 cpp = 1004.70885 Cmaf a verifier.... preff = 101325. pa= 50000. C print*,'PARAMETRES A ENTRER DANS offline.def' print*,'Nombre effectif de traceurs' read(*,*) nq write(*,*) nq print*,'Splitting du pas de temps dynamique' read(*,*) nsplit_dyn write(*,*) nsplit_dyn print*,'Splitting du pas de temps physique' read(*,*) nsplit_phy write(*,*) nsplit_phy print*,'jour0, par rapport au debut du fichier fluxmass' read(*,*) jour0 write(*,*) jour0 print*,'redecoupage (T) ou standard (F)' read(*,'(l1)') redecoupage write(*,*) redecoupage print*,'Lecture fichier start (T)?' read(*,'(l1)') lectstart write(*,'(l1)') lectstart print*,'Traceurs 1 et 2 = Rn222 et Pb' read(*,'(l1)') rnpb write(*,'(l1)') rnpb print*,'Avant T ou arriere F?' read(*,'(l1)') avant write(*,'(l1)') avant print*,'Index min et max des traceurs a initaliser' read(*,*) iqmin,iqmax write(*,*) iqmin,iqmax print*,'localisation (lon, lat, intensite)' if(iqmax.le.0) then iqmin=iqmax+1 endif do iq=iqmin,iqmax read(*,*) lon_stat(iq),lat_stat(iq),intensite(iq) write(*,*) lon_stat(iq),lat_stat(iq),intensite(iq) enddo print*,'Records min et max pour le stokage et couche ' s ,'verticale (specifique reseaux), def 0 0 1' read(*,*) recstkmin,recstkmax,ldec write(*,*) recstkmin,recstkmax,ldec c redefinition du splitting if(mod(nsplit_dyn,nsplit_phy).eq.0) then nsplit=nsplit_phy elseif (mod(nsplit_phy,nsplit_dyn).eq.0) then nsplit=nsplit_dyn else stop'prendre un des deux nsplit doit diviser l autre' endif nsplit_dyn=nsplit_dyn/nsplit nsplit_phy=nsplit_phy/nsplit print*,'nsplit,nsplit_dyn,nsplit_phy' s ,nsplit,nsplit_dyn,nsplit_phy c Lecture eventuelle time_0=0. if(nq.gt.nqmx) stop'surdimensioner nqmx' print*,'av initial0' call initial0(nq*ijp1llm,q) print*,'av lectstart' if(lectstart) then CALL dynetat0("start.nc",nq,vcov,ucov, . teta,q,masse,ps,phis, time_0) print*,'Lecture du start' c On zappe le radon et le plomb. q(:,:,:,1)=q(:,:,:,3) q(:,:,:,2)=q(:,:,:,4) c Initialisation d'un traceur ` 1 pour tester l'impact du lessivage. q(:,:,:,3)=1. else day_ini=0 annee_ref=0 endif day_end=day_ini+nday C print*,'av defrun' c open (99,file='run.def',status='old',form='formatted') clefinit=.true. CALL defrun_new( 99, clefinit, clesphy0 ) c close(99) C print*,'av iniconst' c lecture du jour de demarrage c premiere cnitialisation, eventuellement bidon cALL iniconst CALL inigeom C print*,'ENTREE DANS redecoupenc ou lectfluxnc s pour irec=0' c Lecture des entetes des fichiers de flux if(redecoupage) then call redecoupenc(0,masse,pbaru,pbarv,w,teta,phi, s nrec,avant,airefi,pphis, s t,zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz, s yu1,yv1,ftsol,pctsrf, s frac_impa,frac_nucl,phis) else call lectfluxnc(0,masse,pbaru,pbarv,w,teta,phi, s nrec,avant,airefi,pphis, s t,zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz, s yu1,yv1,ftsol,pctsrf, s frac_impa,frac_nucl,phis) call inigeom endif c controle sur les records if(avant) then recmin=2+(jour0-1)*iday_step recmax=1+(jour0+nday-1)*iday_step else recmin=(jour0-nday)*iday_step+2 recmax=jour0*iday_step+1 endif if(recmin.le.1.or.recmax.ge.nrec) then print*,'Probleme de lecture ' print*,'recmin,recmax,nrec',recmin,recmax,nrec stop endif C on recalcule le nombre de pas de temps dans une journee c istdyn est la frequence de stokage en "pas de temps dynamique" dans c le modele on-line. dtvr est le pas de temps du meme modele en s. c iday-step est le nombre de lectures par jour (grand pas de temps). iday_step=daysec/(istdyn*dtvr) C CALL iniconst c on calcule le pas de temps dtvr = daysec/FLOAT(iday_step) dtphys=dtvr/(nsplit*nsplit_phy) ! a mettre apres iniconst! print*,'Pas de temps physique ',dtphys C call gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi) print *,'Pas de temps dynamique', dtvr print*,'Pas de temps physique ',dtphys c----------------------------------------------------------------------- C maf est ce utile ? call suphec CALL inifilr print*,'Pas de temps physique ',dtphys CALL inidissip( lstardis, nitergdiv, nitergrot, niterh , * tetagdiv, tetagrot , tetatemp ) print*,'Pas de temps physique ',dtphys c----------------------------------------------------------------------- c temps de depart et de fin: c -------------------------- itau = 0 iday = dayref time = time_0 IF(time.GT.1.) THEN time = time-1. iday = iday+1 ENDIF print*,'Pas de temps physique ',dtphys c======================================================================= c Initialisation des fichiers de sortie c======================================================================= print*,'Initialisation de wrgras' do i=1,iip1 print*,'rlonu(',i,')=',rlonu(i) enddo do i=1,jjp1 print*,'rlatu(',i,')=',rlatu(i) enddo do i=1,iip1 print*,'rlonv(',i,')=',rlonv(i) enddo do i=1,jjm print*,'rlatv(',i,')=',rlatv(i) enddo file='sol' call inigrads(1,iip1 s ,rlonv,180./pi,-180.,180.,jjp1,rlatu,-90.,90.,180./pi s ,llm,presnivs,1. s ,dtphys,file,'gcmq2 ') file='atm' print*,'Pas de temps physique ',dtphys call inigrads(2,iip1 s ,rlonv,180./pi,-180.,180.,jjp1,rlatu,-90.,90.,180./pi s ,llm,presnivs,1. s ,dtphys,file,'gcmq2 ') file='pbu' call inigrads(3,iip1 s ,rlonu,180./pi,-180.,180.,jjp1,rlatu,-90.,90.,180./pi s ,llm,presnivs,1. s ,dtphys,file,'gcmq2 ') file='pbv' call inigrads(4,iip1 s ,rlonv,180./pi,-180.,180.,jjm,rlatv,-90.,90.,180./pi s ,llm,presnivs,1. s ,dtphys,file,'gcmq2 ') file='qadv' call inigrads(5,iip1 s ,rlonv,180./pi,-180.,180.,jjp1,rlatu,-90.,90.,180./pi s ,llm,presnivs,1. s ,dtphys,file,'gcmq2 ') print*,'dtphys ap inigrads ',dtphys print*,'Ecriture du fichier de conditions aux limites' open (98,file='limit',form='unformatted') write(98) float(im),float(jm),float(nq),float(nday) write(98) (rlonv(i)*180./pi,i=1,iip1) write(98) (rlatu(j)*180./pi,j=1,jjp1) write(98) ((phis(i,j)/g,i=1,iip1),j=1,jjp1) close(98) CALL dynredem0("restart.nc",day_end,phis,nq) c----------------------------------------------------------------------- c Debut de la boucle en temps: c ---------------------------- do itau=1,nday*iday_step print*,'ITAU ITAU ITAU =',itau injecte='.' stoke='.' c Gestion du temps if (avant) then dtdyn=dtvr/float(nsplit*nsplit_dyn) c irec=1+(jour0-1)*iday_step+itau irec=(jour0-1)*iday_step+itau else dtdyn=-dtvr/float(nsplit*nsplit_dyn) irec=1+jour0*iday_step-itau+1 endif debut=itau.eq.1 c----------------------------------------------------------------------- c Lecture des fichiers fluxmass et physique: c ----------------------------------------------------- print*,'CCCCC Appel a REDECOUPENC ou LECTFLUXNC s au pas itau=',itau if(redecoupage) then call redecoupenc(irec,masse,pbaru,pbarv,w,teta,phi, s nrec,avant,airefi,pphis, s t,zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz, s yu1,yv1,ftsol,pctsrf, s frac_impa,frac_nucl,phis) else call lectfluxnc(irec,masse,pbaru,pbarv,w,teta,phi, s nrec,avant,airefi,pphis, s t,zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz, s yu1,yv1,ftsol,pctsrf, s frac_impa,frac_nucl,phis) endif print*,'TESTPHYS: ON PREND LA PUIS ',1./float(nsplit*nsplit_phy) s ,' DES FRAC A IT=',itau frac_impa(:,:)=frac_impa(:,:)**(1./float(nsplit*nsplit_phy)) frac_nucl(:,:)=frac_nucl(:,:)**(1./float(nsplit*nsplit_phy)) c ... ouverture du fichier de stockage netcdf ... C IF (debut) then dynhistave_file = 'histmoy.nc' day_ini=0 annee_ref=0 t_ops =(1./48.)*daysec t_wrt =(1./48.)*daysec dtav=dtvr/float(nsplit) mode=1 c CALL initdynav(dynhistave_file,day_ini,annee_ref,dtav, c . t_ops, t_wrt, nq,mode, histaveid) pi=2.*asin(1.) do iq=iqmin,iqmax do l=1,llm do j=1,jjp1 do i=1,iip1 q(i,j,l,iq)=0. enddo enddo enddo zlon=lon_stat(iq)*pi/180. zlat=lat_stat(iq)*pi/180. CALL coordij(zlon,zlat,i_stat(iq),j_stat(iq)) enddo ENDIF C calcul de la pression de surface. do j=1,jjp1 do i=1,iip1 smass(i,j)=0. enddo enddo do l=1,llm do j=1,jjp1 do i=1,iip1 smass(i,j)=smass(i,j)+masse(i,j,l) enddo enddo enddo do j=1,jjp1 do i=1,iip1 ps (i,j)=smass(i,j)/aire(i,j)*g end do end do C C calcul de la pression et de pk (fonction d'Exner) C CALL pression ( ip1jmp1, ap, bp, ps, p ) CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) C do jj=1,jjp1 do ii=1,iip1 phis(ii,jj)= phi(ii,jj,1)-teta(ii,jj,1)* s (pks(ii,jj)-pk(ii,jj,1)) end do end do c======================================================================= c TERMES SOURCES OU PUITS do isplit=1,nsplit c======================================================================= c injection par Fred c Print*,Verifier en fct de iapp_tracvl la boucle inj goto 333 c Inj pour le stokage ttes les 0.5h if (itau.gt.32.and.itau.le.56)then print*,'Masse de la maille inj :', s masse(i_stat(1),j_stat(1),1) do iiinj=1,20 q(i_stat(1),j_stat(1),1,1)= s q(i_stat(1),j_stat(1),1,1) s +intensite(1)/masse(i_stat(1),j_stat(1),1)/ s (float(nsplit)*24.*20.) tinj1=tinj1+intensite(1)/(float(nsplit)*24.*20.) enddo print*,'Apres injection a itau ',itau do iq=1,nq write(str2,'(i2.2)') iq call diagadv(q(:,:,:,iq),masse,'Traceur '//str2) enddo print*,'QT TOTALE INJECTEE TR1',tinj1 print*,'fin des injections' endif c Inj pour le stokage ttes les 1h if (itau.gt.16.and.itau.le.28)then print*,'Masse de la maille inj :', s masse(i_stat(1),j_stat(1),1) do iiinj=1,20 q(i_stat(1),j_stat(1),1,1)= s q(i_stat(1),j_stat(1),1,1) s +intensite(1)/masse(i_stat(1),j_stat(1),1)/ s (float(nsplit)*12.*20.) tinj1=tinj1+intensite(1)/(float(nsplit)*12.*20.) enddo print*,'Apres injection a itau ',itau do iq=1,nq write(str2,'(i2.2)') iq call diagadv(q(:,:,:,iq),masse,'Traceur '//str2) enddo print*,'QT TOTALE INJECTEE TR1',tinj1 print*,'fin des injections' endif c Inj pour le stokage ttes les 1.5h if (itau.eq.11.and.isplit.gt.2)then do iiinj=1,20 q(i_stat(1),j_stat(1),1,1)= s q(i_stat(1),j_stat(1),1,1) s +intensite(1)/masse(i_stat(1),j_stat(1),1)/ s (24.*20.) enddo endif if (itau.ge.12.and.itau.le.18)then do iiinj=1,20 q(i_stat(1),j_stat(1),1,1)= s q(i_stat(1),j_stat(1),1,1) s +intensite(1)/masse(i_stat(1),j_stat(1),1)/ s (24.*20.) enddo endif if (itau.eq.19.and.isplit.le.2)then do iiinj=1,20 q(i_stat(1),j_stat(1),1,1)= s q(i_stat(1),j_stat(1),1,1) s +intensite(1)/masse(i_stat(1),j_stat(1),1)/ s (24.*20.) enddo endif c Inj pour le stokage ttes les 4h if (itau.ge.5.and.itau.le.7)then do iiinj=1,20 q(i_stat(1),j_stat(1),1,1)= s q(i_stat(1),j_stat(1),1,1) s +intensite(1)/masse(i_stat(1),j_stat(1),1)/ s (24.*20.) enddo endif c Inj pour le stokage ttes les 3h if (itau.eq.6.and.isplit.gt.2)then do iiinj=1,20 q(i_stat(1),j_stat(1),1,1)= s q(i_stat(1),j_stat(1),1,1) s +intensite(1)/masse(i_stat(1),j_stat(1),1)/ s (24.*20.) enddo endif if (itau.ge.7.and.itau.le.9)then do iiinj=1,20 q(i_stat(1),j_stat(1),1,1)= s q(i_stat(1),j_stat(1),1,1) s +intensite(1)/masse(i_stat(1),j_stat(1),1)/ s (24.*20.) enddo endif if (itau.eq.10.and.isplit.le.2)then do iiinj=1,20 q(i_stat(1),j_stat(1),1,1)= s q(i_stat(1),j_stat(1),1,1) s +intensite(1)/masse(i_stat(1),j_stat(1),1)/ s (24.*20.) enddo endif c Inj pour le stokage ttes les 6h if (itau.eq.3.and.isplit.gt.8)then do iiinj=1,20 q(i_stat(1),j_stat(1),1,1)= s q(i_stat(1),j_stat(1),1,1) s +intensite(1)/masse(i_stat(1),j_stat(1),1)/ s (24.*20.) enddo endif if (itau.eq.4)then do iiinj=1,20 q(i_stat(1),j_stat(1),1,1)= s q(i_stat(1),j_stat(1),1,1) s +intensite(1)/masse(i_stat(1),j_stat(1),1)/ s (24.*20.) enddo endif if (itau.eq.5.and.isplit.le.8)then do iiinj=1,20 q(i_stat(1),j_stat(1),1,1)= s q(i_stat(1),j_stat(1),1,1) s +intensite(1)/masse(i_stat(1),j_stat(1),1)/ s (24.*20.) enddo endif c Inj sur les 4 1ere couches pour le stokage ttes les 3h if (itau.eq.6.and.isplit.gt.2)then do l=1,4 q(i_stat(1),j_stat(1),l,1)= s q(i_stat(1),j_stat(1),l,1) s +intensite(1)/masse(i_stat(1),j_stat(1),l)/ s (24.*4.) enddo endif if (itau.ge.7.and.itau.le.9)then do l=1,4 q(i_stat(1),j_stat(1),l,1)= s q(i_stat(1),j_stat(1),l,1) s +intensite(1)/masse(i_stat(1),j_stat(1),l)/ s (24.*4.) enddo endif if (itau.eq.10.and.isplit.le.2)then do l=1,4 q(i_stat(1),j_stat(1),l,1)= s q(i_stat(1),j_stat(1),l,1) s +intensite(1)/masse(i_stat(1),j_stat(1),l)/ s (24.*4.) enddo endif c Inj pour le stokage ttes les 12h if (itau.eq.2.and.isplit.gt.8)then do iiinj=1,20 q(i_stat(1),j_stat(1),1,1)= s q(i_stat(1),j_stat(1),1,1) s +intensite(1)/masse(i_stat(1),j_stat(1),1)/ s (24.*20.) enddo endif if (itau.eq.3.and.isplit.le.8)then do iiinj=1,20 q(i_stat(1),j_stat(1),1,1)= s q(i_stat(1),j_stat(1),1,1) s +intensite(1)/masse(i_stat(1),j_stat(1),1)/ s (24.*20.) enddo endif c Inj pour le stokage ttes les 24h if (itau.eq.1.and.isplit.gt.32)then q(i_stat(1),j_stat(1),1,1)= s q(i_stat(1),j_stat(1),1,1) s +intensite(1)/masse(i_stat(1),j_stat(1),1)/24. endif if (itau.eq.2.and.isplit.le.8)then q(i_stat(1),j_stat(1),1,1)= s q(i_stat(1),j_stat(1),1,1) s +intensite(1)/masse(i_stat(1),j_stat(1),1)/24. endif c Inj pour le stokage ttes les 3h if (itau.eq.6.and.isplit.gt.2)then do iiinj=1,20 q(i_stat(1),j_stat(1),1,1)= s q(i_stat(1),j_stat(1),1,1) s +intensite(1)/masse(i_stat(1),j_stat(1),1)/ s (24.*20.) enddo endif if (itau.ge.7.and.itau.le.9)then do iiinj=1,20 q(i_stat(1),j_stat(1),1,1)= s q(i_stat(1),j_stat(1),1,1) s +intensite(1)/masse(i_stat(1),j_stat(1),1)/ s (24.*20.) enddo endif if (itau.eq.10.and.isplit.le.2)then do iiinj=1,20 q(i_stat(1),j_stat(1),1,1)= s q(i_stat(1),j_stat(1),1,1) s +intensite(1)/masse(i_stat(1),j_stat(1),1)/ s (24.*20.) enddo endif 333 continue if (itau.eq.6.and.isplit.gt.2)then do l=1,2 q(i_stat(1),j_stat(1),l,1)= s q(i_stat(1),j_stat(1),l,1) s +intensite(1)/masse(i_stat(1),j_stat(1),l)/ s (24.*2.) enddo endif if (itau.ge.7.and.itau.le.9)then do l=1,2 q(i_stat(1),j_stat(1),l,1)= s q(i_stat(1),j_stat(1),l,1) s +intensite(1)/masse(i_stat(1),j_stat(1),l)/ s (24.*2.) enddo endif if (itau.eq.10.and.isplit.le.2)then do l=1,2 q(i_stat(1),j_stat(1),l,1)= s q(i_stat(1),j_stat(1),l,1) s +intensite(1)/masse(i_stat(1),j_stat(1),l)/ s (24.*2.) enddo endif c======================================================================= c FIN DE LA PARTIE TERMES SOURCES c======================================================================= c----------------------------------------------------------------------- c Advection c ---------- print*,'ENTREE DANS VLSPLT' c======================================================================= c do isplit=1,nsplit c======================================================================= CC MAF dans le cas du offline on ne fait pas la difference masse et massem do iii=1,nsplit_dyn nom='pbaru' call wrgrads(3,llm,pbaru(:,:,1),nom,nom) nom='pbarv' call wrgrads(4,llm,pbarv(:,:,1),nom,nom) nom='masse' call wrgrads(3,llm,masse(:,:,1),nom,nom) nom='w' call wrgrads(3,llm,w(:,:,1),nom,nom) print*,'Entree dans vlsplt a itau et au temps',itau,time print*,'Pas advect :',dtdyn print*,'Avant vlsplt' do iq=1,nq write(str2,'(i2.2)') iq call diagadv(q(:,:,:,iq),masse,'Traceur '//str2) enddo do iq=1,nq ctest abder CALL vlsplt(q(1,1,1,iq),2.,masse,w,pbaru,pbarv,dtdyn) CALL vlsplt(q(:,:,:,iq),2.,masse,w,pbaru,pbarv,dtdyn) enddo do iq=iqmin,iqmax write(str2,'(i2.2)') iq nom='q'//str2 call wrgrads(5,1,q(:,:,1,iq),nom,nom) enddo print*,'Apres vlsplt' do iq=1,nq write(str2,'(i2.2)') iq call diagadv(q(:,:,:,iq),masse,'Traceur '//str2) enddo c mise a jour du champ de masse tenant compte de l'advection sur c le pas de temps considere. enddo c----------------------------------------------------------------------- c Physique (lessivage dans phytrac maintenant) MAF c --------------------------------------------------------- if(physic) then if(debut) then latfi(1)=rlatu(1)*180./pi lonfi(1)=0. DO j=2,jjm DO i=1,iim latfi((j-2)*iim+1+i)= rlatu(j)*180./pi lonfi((j-2)*iim+1+i)= rlonv(i)*180./pi ENDDO ENDDO latfi(ngridmx)= rlatu(jjp1)*180./pi lonfi(ngridmx)= 0. endif c C Calcul de paprs et pplay et de temp c ----------------------------------------------------------------- c .... paprs definis aux (llm +1) interfaces des couches .... c .... pplay definis aux ( llm ) milieux des couches .... c ----------------------------------------------------------------- C DO l = 1, llmp1 paprs( 1,l ) = p(1,1,l) ig0 = 2 DO j = 2, jjm DO i =1, iim paprs( ig0,l ) = p(i,j,l) ig0 = ig0 +1 ENDDO ENDDO paprs( ngridmx,l ) = p(1,jjp1,l) ENDDO c unskap = 1./ kappa DO l=1,llm pksurcp = pk(1,1,l) / cpp pplay(1,l) = preff * pksurcp ** unskap tempfi(1,l) = teta(1,1,l) * pksurcp ig0 = 2 DO j = 2, jjm DO i = 1, iim pksurcp = pk(i,j,l) / cpp pplay(ig0,l) = preff * pksurcp ** unskap tempfi(ig0,l) = teta(i,j,l) * pksurcp ig0 = ig0 + 1 ENDDO ENDDO pksurcp = pk(1,jjp1,l) / cpp pplay(ig0,l) = preff * pksurcp ** unskap tempfi (ig0,l) = teta(1,jjp1,l) * pksurcp ENDDO c On passe le traceur et le geopot de surf sur la grille physique call gr_dyn_fi(llm*nq,iip1,jjp1,ngridmx,q,qfi) c call gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,pphis) do iii=1,nsplit_phy C lafin=.false. ! en attendant mieux. print*,'dtphys avant phytrac ',dtphys print*,'TESTPHYS: APPEL A PHYTRAC IT=',itau call phytrac(rnpb, I debutphy, lafin, I nq, I ngridmx,llm,dtphys, I t,paprs,pplay, I zmfu, zmfd, zen_u, zde_u, zen_d, zde_d, I coefkz,yu1,yv1,ftsol,pctsrf,latfi, I frac_impa, frac_nucl, I lonfi,presnivs,airefi,pphis, O qfi) C debutphy= .FALSE. C enddo ! nsplit_phy C c On passe le traceur sur la grille dynamique. call gr_fi_dyn(llm*nq,ngridmx,iip1,jjp1,qfi,q) endif itauav=(itau-1)*nsplit+isplit c itauav=itau*nsplit+isplit c CALL writedynav(histaveid, nq,mode, itauav,vcov , c , ucov,teta,pk,phi,q,masse,ps,phis) c qmoy(:,:,:)=qmoy(:,:,:)+q(:,:,1,:) do iq=1,nq do j=1,jjp1 do i=1,iip1 qmoy(i,j,iq)=qmoy(i,j,iq)+q(i,j,1,iq) enddo enddo enddo enddo ! isplit=1,nsplit call histsync c----------------------------------------------------------------------- c Calcul de ucov et vcov pour les sorties c ------------------------------------- call massbar(masse, massebx, masseby ) do l=1,llm do j=1,jjp1 do i=1,iip1 ucov(i,j,l)=pbaru(i,j,l)/massebx(i,j,l)*cu(i,j)*cu(i,j) s /istdyn enddo enddo do j=1,jjm do i=1,iip1 vcov(i,j,l)=pbarv(i,j,l)/masseby(i,j,l)*cv(i,j)*cv(i,j) s /istdyn enddo enddo enddo iday= dayref+itau/iday_step time= s float(itau-(iday-dayref)*iday_step)/iday_step+time_0 IF(time.GT.1.) THEN time = time-1. iday = iday+1 ENDIF C c======================================================================= do iq=iqmin,iqmax write(str2,'(i2.2)') iq nom='q'//str2 call wrgrads(1,1,qmoy(:,:,iq),nom,nom) call wrgrads(2,1,q(:,:,ldec,iq),nom,nom) enddo do iq=1,nq do j=1,jjp1 do i=1,iip1 qmoy(i,j,iq)=0. enddo enddo enddo print*,'Record=',irec,' stoke=',stoke, s ' injecte=',injecte enddo ! itau CALL dynredem1("restart.nc",0.0, . vcov,ucov,teta,q,nq,masse,ps) c======================================================================= C fermeture des fichiers netcdf CALL histclo C 300 FORMAT('1',/,15x,'run du pas',i7,2x,'au pas',i7,2x, s 'c"est a dire du jour',i7,3x,'au jour',i7,/,/) C---------------------------------------------------------------------- stop end