PROGRAM offline USE ioipsl IMPLICIT NONE c ...... Version du 17/04/96 .......... c======================================================================= c c Auteur: P. Le Van /L. Fairhead/F.Hourdin c ------- c Modif special traceur F.Forget 05/94 c c Modifs M-A Filiberti 07/99 C c Objet: c ------ c c GCM LMD nouvelle grille c c======================================================================= c ... modification de l'integration de q ( 26/04/94 ) .... c----------------------------------------------------------------------- c Declarations: c ------------- #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 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 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 integer jour0,isplit,nsplit_dyn,nsplit_phy,nsplit logical debut,lectstart,rnpb 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) c----------------------------------------------------------------------- c Initialisations: c ---------------- descript = 'Run GCM LMDZ' c----------------------------------------------------------------------- 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*,'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*,'Injection sur combien de pas de temps' read(*,*) npasinject write(*,*) npasinject 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' else day_ini=0 anne_ini=0 endif day_end=day_ini+nday print*,'av defrun' c open (99,file='run.def',status='old',form='formatted') clefinit=.true. CALL defrun_new( 99, clefinit, clesphy0 ) c close(99) print*,'av iniconst' c lecture du jour de demarrage c premiere initialisation, eventuellement bidon CALL iniconst CALL inigeom C c Lecture des entetes des fichiers de flux if(redecoupage) then call redecoupe(0,masse,pbaru,pbarv,w,teta,phi, 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) else call lectflux(0,masse,pbaru,pbarv,w,teta,phi, 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,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' file='sol' call inigrads(1,iip1 s ,rlonv,180./pi,-190.,190.,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,-190.,190.,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,anne_ini,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) irec=1+(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 redecoupe au pas itau=',itau if(redecoupage) then call redecoupe(irec,masse,pbaru,pbarv,w,teta,phi, 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) else call lectflux(irec,masse,pbaru,pbarv,w,teta,phi, 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,phis) endif c ... ouverture du fichier de stockage netcdf ... C IF (debut) then dynhistave_file = 'dyn_hist_ave.nc' day_ini=0 anne_ini=0 t_ops = periodav * daysec t_wrt = periodav * daysec dtav=dtvr/float(nsplit) mode=1 CALL initdynav(dynhistave_file,day_ini,anne_ini,dtav, . t_ops, t_wrt, nq, 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 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 c======================================================================= c terme source reparti sur npasinject pas c On injecte tous les traceurs de la meme facon if (itau.le.npasinject) then injecte='X' do iq=iqmin,iqmax print*,'INJECTION AU POINT ',lon_stat(iq),lat_stat(iq) s ,' I,J=',i_stat(iq),j_stat(iq),jjp1-j_stat(iq)+1 q(i_stat(iq),j_stat(iq),1,iq)= s q(i_stat(iq),j_stat(iq),1,iq) s +intensite(iq)/masse(i_stat(iq),j_stat(iq),1)/npasinject enddo print*,'fin des injections' endif c======================================================================= c FIN DE LA PARTIE TERMES SOURCES c======================================================================= c----------------------------------------------------------------------- c Advection 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 do iq=1,nq CALL vlsplt(q(:,:,:,iq),2.,masse,w,pbaru,pbarv,dtdyn) enddo c mise à jour du champd de masse tenant compte de l'advection sur c le pas de temps considere. enddo do iq=1,nq write(str2,'(i2.2)') iq if(mod(itau,1).eq.0) s call diagadv(q(:,:,:,iq),masse,'Traceur '//str2) 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) call gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,pphis) do iii=1,nsplit_phy C print*,'dtphys avant phytrac ',dtphys call phytrac(rnpb, I ecritphy, debutphy, I nq, I ngridmx,llm,dtphys, I tempfi,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 CALL writedynav(histaveid, nq, itauav,vcov , , ucov,teta,pk,phi,q,masse,ps,phis) 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======================================================================= if(irec.ge.recstkmin.and.irec.le.recstkmax) then if (mod(itau,4).eq.0) then stoke='X' do iq=1,nq write(str2,'(i2.2)') iq nom='q'//str2 call wrgrads(1,1,q(:,:,:,iq),nom,nom) call wrgrads(2,1,q(:,:,ldec,iq),nom,nom) enddo endif endif 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 dujour',i7,3x,'au jour',i7,/,/) end C----------------------------------------------------------------------