Index: DZ.3.3/branches/rel-LF/libf/dyn3d/fluxstoke.F =================================================================== --- /LMDZ.3.3/branches/rel-LF/libf/dyn3d/fluxstoke.F (revision 224) +++ (revision ) @@ -1,149 +1,0 @@ - SUBROUTINE fluxstoke(pbaru,pbarv,masse,teta,phi,phis) -c -c Auteur : F. Hourdin -c -c -ccc .. Modif. P. Le Van ( 20/12/97 ) ... -c - IMPLICIT NONE -c -#include "dimensions.h" -#include "paramet.h" -#include "comconst.h" -#include "comvert.h" -#include "comgeom.h" -#include "tracstoke.h" - - - REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) - REAL masse(ip1jmp1,llm),teta(ip1jmp1,llm),phi(ip1jmp1,llm) - REAL phis(ip1jmp1) - - REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm) - REAL massem(ip1jmp1,llm),tetac(ip1jmp1,llm),phic(ip1jmp1,llm) - - REAL pbarug(ip1jmp1,llm),pbarvg(iip1,jjm,llm),wg(ip1jmp1,llm) - - REAL pbarvst(iip1,jjp1,llm) - - - INTEGER iadvtr - INTEGER ij,l,irec,i,j - - SAVE iadvtr, massem,pbaruc,pbarvc,irec - SAVE phic,tetac - logical first - save first - data first/.true./ - DATA iadvtr/0/ - - if(first) then -#ifdef CRAY - CALL ASSIGN("assign -N ieee -F null f:fluxmass") -#endif - open(47,file='fluxmass',form='unformatted', - s access='direct',recl=4*(6*ijp1llm)) - irec=1 - first=.false. - - open(77,file='fluxmass.ctl',status='unknown',form='formatted') - - endif - - - IF(iadvtr.EQ.0) THEN - CALL initial0(ijp1llm,phic) - CALL initial0(ijp1llm,tetac) - CALL initial0(ijp1llm,pbaruc) - CALL initial0(ijmllm,pbarvc) - ENDIF - -c accumulation des flux de masse horizontaux - DO l=1,llm - DO ij = 1,ip1jmp1 - pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l) - tetac(ij,l) = tetac(ij,l) + teta(ij,l) - phic(ij,l) = phic(ij,l) + phi(ij,l) - ENDDO - DO ij = 1,ip1jm - pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l) - ENDDO - ENDDO - -c selection de la masse instantannee des mailles avant le transport. - IF(iadvtr.EQ.0) THEN - CALL SCOPY(ip1jmp1*llm,masse,1,massem,1) - ENDIF - - iadvtr = iadvtr+1 - - -c Test pour savoir si on advecte a ce pas de temps - IF ( iadvtr.EQ.istdyn ) THEN - -c normalisation - DO l=1,llm - DO ij = 1,ip1jmp1 - pbaruc(ij,l) = pbaruc(ij,l)/float(istdyn) - tetac(ij,l) = tetac(ij,l)/float(istdyn) - phic(ij,l) = phic(ij,l)/float(istdyn) - ENDDO - DO ij = 1,ip1jm - pbarvc(ij,l) = pbarvc(ij,l)/float(istdyn) - ENDDO - ENDDO - -c traitement des flux de masse avant advection. -c 1. calcul de w -c 2. groupement des mailles pres du pole. - - CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg ) - - do l=1,llm - do j=1,jjm - do i=1,iip1 - pbarvst(i,j,l)=pbarvg(i,j,l) - enddo - enddo - do i=1,iip1 - pbarvst(i,jjp1,l)=0. - enddo - enddo - - iadvtr=0 - - irec=irec+1 - write(47,rec=1) float(irec),dtvr,float(istdyn), - s float(iim),float(jjm),float(llm),rlonu,rlonv,rlatu,rlatv - s ,aire,phis - write(47,rec=irec) massem,pbarug,pbarvst,wg,tetac,phic - -c on reinitialise a zero les flux de masse cumules. - - write(77,'(a4,2x,a40)') - & 'DSET ','^fluxmass' - - write(77,'(a12)') 'UNDEF 1.0E30' - write(77,'(a5,1x,a40)') 'TITLE ','Titre a voir' - call formcoord(77,iip1,rlonv,180./pi,.false.,'XDEF') - call formcoord(77,jjp1,rlatu,180./pi,.true.,'YDEF') - call formcoord(77,llm,presnivs,1.,.false.,'ZDEF') - write(77,'(a4,i10,a30)') - & 'TDEF ',irec,' LINEAR 02JAN1987 1DY ' - write(77,'(a4,2x,i5)') 'VARS',6 - write(77,1000) 'masse',llm,99,'masse ' - write(77,1000) 'pbaru',llm,99,'pbaru ' - write(77,1000) 'pbarv',llm,99,'pbarv ' - write(77,1000) 'w ',llm,99,'w ' - write(77,1000) 'teta ',llm,99,'teta ' - write(77,1000) 'phi ',llm,99,'phi ' - write(77,'(a7)') 'ENDVARS' - -1000 format(a5,3x,i4,i3,1x,a39) - - - - ENDIF ! if iadvtr.EQ.istdyn - - RETURN - END Index: DZ.3.3/branches/rel-LF/libf/dyn3d/lectflux.F =================================================================== --- /LMDZ.3.3/branches/rel-LF/libf/dyn3d/lectflux.F (revision 224) +++ (revision ) @@ -1,233 +1,0 @@ - SUBROUTINE lectflux(irec,massem,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,phis) - - 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 ngridmx,ig,l - parameter (ngridmx=iim*(jjm-1)+2) - INTEGER nbsrf - PARAMETER (nbsrf=4) ! nombre de sous-fractions pour une maille - - real zmfd(ngridmx,llm),zde_d(ngridmx,llm),zen_d(ngridmx,llm) - 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 yu1(ngridmx), yv1(ngridmx) - real ftsol(ngridmx,nbsrf),pctsrf(ngridmx,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(ngridmx,llm) - - real massem(ip1jmp1,llm),tetan(ip1jmp1,llm) - real pbarun(iip1,jjp1,llm),pbarvn(iip1,jjm,llm) - real pbarvst(iip1,jjp1,llm) - real wn(iip1,jjp1,llm),phin(iip1,jjp1,llm) - real phis(iip1,jjp1) - - real airefi(ngridmx) - - real xlecn(ngridmx,9*llm+2+2*nbsrf) - - real zcontrole(ngridmx),zmass,tmpdyn(iip1,jjp1),zflux - - real ziadvtrac,zrec,ziadvtrac2,zrec2 - real zim,zjm,zlm,zklon,zklev - - real zpi - - 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 ouverture des fichiers -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*ijp1llm)) - read(47,rec=1) zrec,dtvr,ziadvtrac,zim,zjm,zlm, - s rlonu,rlonv,rlatu,rlatv,aire - s ,phis - print*,'zrec,dtvr,ziadvtrac,zim,zjm,zlm' - print*,zrec,dtvr,ziadvtrac,zim,zjm,zlm - print*,rlonv - - - -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*ngridmx*(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 - print*,'nrec=',nrec - - istdyn=ziadvtrac - istphy=ziadvtrac2 - -c================================================================== -c Fin des initialisations - else ! irec=0 -c================================================================== - - -c----------------------------------------------------------------------- -c Lecture des fichiers fluxmass et physique: -c ----------------------------------------------------- - -c Variables dynamiques - read(47,rec=irec) massem,pbarun,pbarvst,wn,tetan,phin - do l=1,llm - do j=1,jjm - do i=1,iip1 - pbarvn(i,j,l)=pbarvst(i,j,l) - enddo - enddo - enddo - -c Variables physiques - read(49,rec=irec) ((xlecn(ig,l),ig=1,ngridmx), - s l=1,9*llm+2+2*nbsrf) - - do l=1,llm - do ig=1,ngridmx - 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,ngridmx - ftsol(ig,l)=xlecn(ig,itsol+l-1) - pctsrf(ig,l)=xlecn(ig,ipsf+l-1) - enddo - enddo - do ig=1,ngridmx - 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,ngridmx - 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,ngridmx - 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----------------------------------------------------------------------- - - print*,'Ap redec irec' - - call gr_dyn_fi(llm,iip1,jjp1,ngridmx,massem,massefi) - - do ig=1,ngridmx - zcontrole(ig)=1. - enddo -c zmass=(max(massem(ig,l),massem(ig,l-1))/airefi(ig) - do l=2,llm - do ig=1,ngridmx - 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,ngridmx - 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,ngridmx,iip1,jjp1,zcontrole,tmpdyn) - - do l=1,llm - do ig=1,ngridmx - 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 - - Index: DZ.3.3/branches/rel-LF/libf/dyn3d/offline.F =================================================================== --- /LMDZ.3.3/branches/rel-LF/libf/dyn3d/offline.F (revision 224) +++ (revision ) @@ -1,640 +1,0 @@ - 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---------------------------------------------------------------------- Index: DZ.3.3/branches/rel-LF/libf/dyn3d/redecoupe.F =================================================================== --- /LMDZ.3.3/branches/rel-LF/libf/dyn3d/redecoupe.F (revision 224) +++ (revision ) @@ -1,637 +1,0 @@ - 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 - - Index: DZ.3.3/branches/rel-LF/libf/dyn3d/run.def =================================================================== --- /LMDZ.3.3/branches/rel-LF/libf/dyn3d/run.def (revision 224) +++ (revision ) @@ -1,187 +1,0 @@ - ------------------------------------------------------------------------ -Parametres de controle du run: ------------------------------- - -- Jour de l'etat initial ( = 350 si 20 Decembre ,par expl. ,comme ici ) - dayref - 350 - -- Annee de l'etat initial ( avec 4 chiffres ) - anneeref - 1987 - -- Nombre de jours d'integration - nday - 2 - -- nombre de pas par jour (multiple de iperiod) ( ici pour dt = 1 min ) - day_step - 1440 - -- periode pour le pas Matsuno (en pas) - iperiod - 5 - -- periode de sortie des variables de controle (en pas) - iconser - 5 - -- periode d'ecriture du fichier histoire (en jour) - iecri - 100 - -- periode de stockage fichier histmoy (en jour) - periodav - 60. - -- periode de la dissipation (en pas) - idissip - 10 - -- choix de l'operateur de dissipation (star ou non star ) - lstardis - T - -- nombre d'iterations de l'operateur de dissipation gradiv -nitergdiv - 1 - -- nombre d'iterations de l'operateur de dissipation nxgradrot -nitergrot - 2 - -- nombre d'iterations de l'operateur de dissipation divgrad - niterh - 2 - -- temps de dissipation des plus petites long.d ondes pour u,v (gradiv) - tetagdiv - 3000. - -- temps de dissipation des plus petites long.d ondes pour u,v(nxgradrot) - tetagrot - 3000. - -- temps de dissipation des plus petites long.d ondes pour h ( divgrad) - tetatemp - 3000. - -- coefficient pour gamdissip - coefdis - 0. - -- choix du shema d'integration temporelle (Matsuno ou Matsuno-leapfrog) - purmats - F - -- avec ou sans physique - physic - T - -- periode de la physique (en pas) - iphysiq - 15 - -- frequence (en jours ) de l'ecriture du fichier histphy - ecritphy - 4 - -- Cycle diurne ou non - cycle_diurne - T - -- Soil Model ou non - soil_model - F - -- Choix ou non de New oliq - new_oliq - T - -- Orodr ou non pour l orographie - ok_orodr - T - -- Orolf ou non pour l orographie - ok_orolf - T - -- Si = .T. , lecture du fichier limit avec la bonne annee - ok_limitvrai - T - -- Nombre d'appels des routines de rayonnements ( par jour) - nbapp_rad - 12 - -- Flag pour la convection ( 1 pour LMD, 2 pour Tiedtke, 3 CCM ) - iflag_con - 2 - -- longitude en degres du centre du zoom - clon - 63. - -- latitude en degres du centre du zoom - clat - 0. - -- facteur de grossissement du zoom,selon longitude - grossismx - 2.5 - -- facteur de grossissement du zoom ,selon latitude - grossismy - 2.5 - -- Fonction f(y) hyperbolique si = .true. , sinon sinusoidale - fxyhypb - T - -- extension en longitude de la zone du zoom ( fraction de la zone totale) - dzoomx - 0.291667 - -- extension en latitude de la zone du zoom ( fraction de la zone totale) - dzoomy - 0.291667 - -- Fonction f(y) avec y = Sin(latit.) si = .true. , sinon y = latit. - ysinus - T - -c -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -cc -c Autres parametres ou cles , non places dans run.def mais -c dans les routines . -cc -ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -cc -cc A ) iadv(iq ) , iq =1, nqmx , avec iadv = 1 , shema humid. -cc specifique LMD -cc -cc avec iadv = 2 , shema amont -cc -cc ( avec a l'interieur de amont , la cle amont_qs ) -cc -cc avec iadv = 3 , shema Van_Leer -cc -cc dans *** gcm.F **** -cc -cc B ) ok_oasis , pour coupleur Oasis , dans *** physiq.F *** -cc -cc C ) ok_ocean dans *** physiq.F *** -cc -cc D ) ok_journe , ok_instan, ok_region pour sorties analyses -c -cc dans *** physiq.F *** -cc -cc E ) nvm : nbre de vegetations -cc dans *** physiq.F *** - -cc F ) deltay : deplacement ( en degres ) en y de la zone du zoom -cc dans *** inigeom.F *** -