! c c SUBROUTINE phytrac (rnpb, I nstep, I julien, I gmtime, I debutphy, I lafin, I nqmax, I nlon, I nlev, I pdtphys, I u, I v, I t_seri, I paprs, I pplay, I pmfu, I pmfd, I pen_u, I pde_u, I pen_d, I pde_d, I cdragh, I coefh, I fm_therm, I entr_therm, I yu1, I yv1, I ftsol, I pctsrf, I xlat, I frac_impa, I frac_nucl, I xlon, I presnivs, I pphis, I pphi, I albsol, I sh, I rh, I cldfra, I rneb, I diafra, I cldliq, I itop_con, I ibas_con, I pmflxr, I pmflxs, I prfl, I psfl, I da, I phi, I mp, I upwd, I dnwd, I aerosol_couple, I flxmass_w, I tau_inca, I piz_inca, I cg_inca, I ccm, I rfname, O tr_seri) USE ioipsl USE dimphy USE mod_grid_phy_lmdz USE mod_phys_lmdz_para USE comgeomphy USE iophy USE vampir IMPLICIT none c====================================================================== c Auteur(s) FH c Objet: Moniteur general des tendances traceurs c cAA Remarques en vrac: cAA-------------------- cAA 1/ le call phytrac se fait avec nqmax-2 donc nous avons bien cAA les vrais traceurs (nbtr) dans phytrac (pas la vapeur ni eau liquide) cAA 2/ Le choix du radon et du pb se fait juste avec un data cAA (peu propre). Peut-etre pourrait-on prevoir dans l'avenir cAA une variable "type de traceur" c====================================================================== #include "YOMCST.h" #include "dimensions.h" #include "indicesol.h" #include "clesphys.h" #include "temps.h" #include "paramet.h" #include "control.h" #include "advtrac.h" #include "thermcell.h" c====================================================================== c Arguments: c c EN ENTREE: c ========== c c divers: c ------- c integer nlon ! nombre de points horizontaux integer nlev ! nombre de couches verticales integer nqmax ! nombre de traceurs auxquels on applique la physique integer nstep ! appel physique integer julien !jour julien integer itop_con(nlon) integer ibas_con(nlon) real gmtime real pdtphys ! pas d'integration pour la physique (seconde) real t_seri(nlon,nlev) ! temperature real tr_seri(nlon,nlev,nbtr) ! traceur real u(klon,klev) real v(klon,klev) real sh(nlon,nlev) ! humidite specifique real rh(nlon,nlev) ! humidite relative real cldliq(nlon,nlev) ! eau liquide nuageuse real cldfra(nlon,nlev) ! fraction nuageuse (tous les nuages) real diafra(nlon,nlev) ! fraction nuageuse (convection ou stratus artificiels) real rneb(nlon,nlev) ! fraction nuageuse (grande echelle) real albsol(nlon) ! albedo surface real paprs(nlon,nlev+1) ! pression pour chaque inter-couche (en Pa) real ps(nlon) ! pression surface real pplay(nlon,nlev) ! pression pour le mileu de chaque couche (en Pa) real pphi(nlon,klev) ! geopotentiel real pphis(klon) REAL presnivs(klev) logical debutphy ! le flag de l'initialisation de la physique logical lafin ! le flag de la fin de la physique c Olivia integer nsplit REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1) !--lessivage convection REAL prfl(klon,klev+1), psfl(klon,klev+1) !--lessivage large-scale LOGICAL aerosol_couple REAL flxmass_w(klon,klev) CHARACTER(len=8) :: solsym(nqmax) integer la REAL :: tau_inca(klon,klev,9,2) REAL :: piz_inca(klon,klev,9,2) REAL :: cg_inca(klon,klev,9,2) character*4 :: rfname(9) REAL :: ccm(klon,klev,2) c c convection: c ----------- c REAL pmfu(nlon,nlev) ! flux de masse dans le panache montant REAL pmfd(nlon,nlev) ! flux de masse dans le panache descendant REAL pen_u(nlon,nlev) ! flux entraine dans le panache montant c c thermiques: c ----------- c real fm_therm(klon,klev+1),entr_therm(klon,klev) real fm_therm1(klon,klev) c REAL pde_u(nlon,nlev) ! flux detraine dans le panache montant REAL pen_d(nlon,nlev) ! flux entraine dans le panache descendant REAL pde_d(nlon,nlev) ! flux detraine dans le panache descendant c KE real da(nlon,nlev),phi(nlon,nlev,nlev),mp(nlon,nlev) REAL upwd(nlon,nlev) ! saturated updraft mass flux REAL dnwd(nlon,nlev) ! saturated downdraft mass flux c c Couche limite: c -------------- c REAL cdragh(nlon,nlev)! coeff drag pour T et Q REAL coefh(nlon,nlev) ! coeff melange CL REAL yu1(nlon) ! vents au premier niveau REAL yv1(nlon) ! vents au premier niveau REAL xlat(nlon) ! latitudes pour chaque point REAL xlon(nlon) ! longitudes pour chaque point c c Lessivage: c ---------- c c pour le ON-LINE c REAL frac_impa(nlon,nlev) ! fraction d'aerosols impactes REAL frac_nucl(nlon,nlev) ! fraction d'aerosols nuclees c cAA cAA Arguments necessaires pour les sources et puits de traceur: cAA ---------------- cAA real ftsol(nlon,nbsrf) ! Temperature du sol (surf)(Kelvin) real pctsrf(nlon,nbsrf) ! Pourcentage de sol f(nature du sol) c abder real pftsol1(nlon),pftsol2(nlon),pftsol3(nlon),pftsol4(nlon) real ppsrf1(nlon),ppsrf2(nlon),ppsrf3(nlon),ppsrf4(nlon) c fin cAA ---------------------------- cAA VARIABLES LOCALES TRACEURS cAA ---------------------------- cAA cAA Sources et puits des traceurs: cAA ------------------------------ cAA cAA Pour l'instant seuls les cas du rn et du pb ont ete envisages. REAL source(klon,nqmax) ! a voir lorsque le flux est prescrit cAA cAA Pour la source de radon et son reservoir de sol cAA ................................................ REAL,save,allocatable :: trs(:,:) ! Conc. radon ds le sol c$OMP THREADPRIVATE(trs) cym SAVE trs REAL :: trs_tmp(klon_glo) REAL,save,allocatable :: masktr(:,:) ! Masque reservoir de sol traceur c Masque de l'echange avec la surface c (1 = reservoir) ou (possible => 1 ) c$OMP THREADPRIVATE(masktr) cym SAVE masktr REAL,save,allocatable :: fshtr(:,:) ! Flux surfacique dans le reservoir de sol c$OMP THREADPRIVATE(fshtr) cym SAVE fshtr REAL,save,allocatable :: hsoltr(:) ! Epaisseur equivalente du reservoir de sol c$OMP THREADPRIVATE(hsoltr) cym SAVE hsoltr REAL,save,allocatable :: tautr(:) ! Constante de decroissance radioactive c$OMP THREADPRIVATE(tautr) cym SAVE tautr REAL,save,allocatable :: vdeptr(:) ! Vitesse de depot sec dans la couche Brownienne c$OMP THREADPRIVATE(vdeptr) cym SAVE vdeptr REAL,save,allocatable :: scavtr(:) ! Coefficient de lessivage c$OMP THREADPRIVATE(scavtr) cym SAVE scavtr cAA CHARACTER*2 itn C maf ioipsl CHARACTER*2 str2 INTEGER nhori, nvert REAL zsto, zout, zjulian INTEGER, SAVE :: nid_tra c$OMP THREADPRIVATE(nid_tra) INTEGER, SAVE :: nid_tra2,nid_tra3 c$OMP THREADPRIVATE(nid_tra2,nid_tra3) INTEGER ndex(1) INTEGER ndex2d(iim*(jjm+1)),ndex3d(iim*(jjm+1)*klev) REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D REAL zx_tmp_2d(iim,jjm+1), zx_tmp_3d(iim,jjm+1,klev) REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1) c integer itau_w ! pas de temps ecriture = nstep + itau_phy c logical ok_sync parameter (ok_sync = .true.) C C nature du traceur c logical,save,allocatable :: aerosol(:) ! Nature du traceur c ! aerosol(it) = true => aerosol c ! aerosol(it) = false => gaz c ! nat_trac(it) = 1. aerosol logical,save,allocatable :: clsol(:) ! clsol(it) = true => CL sol calculee logical,save,allocatable :: radio(:) ! radio(it)=true => decroisssance radioactive c$OMP THREADPRIVATE(aerosol,clsol,radio) cym save aerosol,clsol,radio C c====================================================================== c c Declaration des procedures appelees c c--modif convection tiedtke INTEGER i, k, it INTEGER iq, iiq REAL delp(klon,klev) c--end modif c c Variables liees a l'ecriture de la bande histoire physique c c Variables locales pour effectuer les appels en serie c---------------------------------------------------- c REAL d_tr(klon,klev), d_trs(klon) ! tendances de traceurs REAL d_tr_cl(klon,klev,nbtr) ! tendance de traceurs couche limite REAL d_tr_cv(klon,klev,nbtr) ! tendance de traceurs conv pour chq traceur REAL d_tr_th(klon,klev,nbtr) ! la tendance des thermiques REAL d_tr_dec(klon,klev,nbtr) ! la tendance de la decroissance c ! radioactive du rn - > pb REAL d_tr_lessi_impa(klon,klev,nbtr) ! la tendance du lessivage c ! par impaction REAL d_tr_lessi_nucl(klon,klev,nbtr) ! la tendance du lessivage c ! par nucleation REAL fluxrn(klon,klev) REAL fluxpb(klon,klev) REAL pbimpa(klon,klev) REAL pbnucl(klon,klev) REAL rn(klon,klev) REAL pb(klon,klev) REAL flestottr(klon,klev,nbtr) ! flux de lessivage c ! dans chaque couche real zmasse(klon,klev) real ztra_th(klon,klev) C character*20 modname character*80 abort_message c c Controles c------------- logical first,couchelimite,convection,lessivage,sorties, s rnpb,inirnpb save first,couchelimite,convection,lessivage, s sorties,inirnpb c$OMP THREADPRIVATE(first,couchelimite,convection,lessivage, c$OMP+ sorties,inirnpb) c data first,couchelimite,convection,lessivage,sorties c s /.true.,.true.,.false.,.true.,.true./ c Olivia data first,couchelimite,convection,lessivage, s sorties s /.true.,.true.,.true.,.true.,.true./ ! Variables needed for configuration with INCA INTEGER :: lastgas INTEGER :: ncsec REAL, PARAMETER :: dry_mass = 28.966 REAL, POINTER :: hbuf(:) REAL, ALLOCATABLE :: obuf(:) REAL :: calday REAL :: pdel(klon,klev) c c====================================================================== modname='phytrac' ps(:)=paprs(:,1) if (debutphy) then allocate( trs(klon,nbtr) ) allocate( masktr(klon,nbtr)) allocate( fshtr(klon,nbtr) ) allocate( hsoltr(nbtr)) allocate( tautr(nbtr)) allocate( vdeptr(nbtr)) allocate( scavtr(nbtr)) allocate( aerosol(nbtr)) allocate( clsol(nbtr)) allocate( radio(nbtr)) ! FH 2008/05/09 correction de la frequence d'ecriture des traceurs ! ecrit_tra = FLOAT(NINT(86400./pdtphys *ecritphy)) print*,'dans phytrac ',pdtphys,ecrit_tra if(nbtr.lt.nqmax) then c print*,'NQMAX=',nqmax c print*,'NBTR=',nbtr abort_message='See above' call abort_gcm(modname,abort_message,1) endif inirnpb=rnpb PRINT*, 'La frequence de sortie traceurs est ', ecrit_tra C c============================================================= c Initialisation des sorties c============================================================= #ifdef CPP_IOIPSL #include "ini_histrac.h" #endif c====================================================================== c Initialisation de certaines variables pour le Rn et le Pb c====================================================================== c Initialisation du traceur dans le sol (couche limite radonique) c c print*,'valeur de debut dans phytrac :',debutphy trs(:,:) = 0. c$OMP MASTER if (is_mpi_root) then trs_tmp(:)=0. open (99,file='starttrac',status='old', . err=999,form='formatted') read(99,*) (trs_tmp(i),i=1,klon_glo) 999 close(99) endif c$OMP END MASTER call Scatter(trs_tmp,trs(:,1)) c print*, 'apres starttrac' c Initialisation de la fraction d'aerosols lessivee c d_tr_lessi_impa(:,:,:) = 0. d_tr_lessi_nucl(:,:,:) = 0. c c Initialisation de la nature des traceurs c DO it = 1, nqmax aerosol(it) = .FALSE. ! Tous les traceurs sont des gaz par defaut radio(it) = .FALSE. ! Par defaut pas de passage par radiornpb clsol(it) = .FALSE. ! Par defaut couche limite avec flux prescrit ENDDO c ENDIF ! fin debutphy c Initialisation du traceur dans le sol (couche limite radonique) if(inirnpb) THEN c radio(1)= .true. radio(2)= .true. clsol(1)= .true. clsol(2)= .true. aerosol(2) = .TRUE. ! le Pb est un aerosol c call initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr . ,vdeptr,scavtr) inirnpb=.false. endif IF (config_inca == 'none') THEN DO i=1,nlon pftsol1(i) = ftsol(i,1) pftsol2(i) = ftsol(i,2) pftsol3(i) = ftsol(i,3) pftsol4(i) = ftsol(i,4) ppsrf1(i) = pctsrf(i,1) ppsrf2(i) = pctsrf(i,2) ppsrf3(i) = pctsrf(i,3) ppsrf4(i) = pctsrf(i,4) ENDDO ELSE ! config_inca /=none #ifdef INCA call VTe(VTphysiq) call VTb(VTinca) !====================================================================== ! Chimie !====================================================================== calday = FLOAT(julien) + gmtime ncsec = NINT (86400.*gmtime) DO k = 1, nlev pdel(:,k) = paprs(:,k) - paprs (:,k+1) END DO IF (config_inca == 'aero') THEN CALL aerosolmain (aerosol_couple, $ tr_seri, $ pdtphys, $ pplay, $ pdel, $ prfl, $ pmflxr, $ psfl, $ pmflxs, $ pmfu, $ itop_con, $ ibas_con, $ pphi, $ airephy, ! paire, $ nstep, $ rneb, ! for chimiaq $ t_seri, ! for chimiaq $ rh, $ tau_inca, $ piz_inca, $ cg_inca, $ rfname, $ ccm, $ lafin) END IF CALL chemmain (tr_seri, !mmr $ nstep, !nstep $ calday, !calday $ julien, !ncdate $ ncsec, !ncsec $ 1, !lat $ pdtphys, !delt $ paprs(1,1), !ps $ pplay, !pmid $ pdel, !pdel $ airephy, $ pctsrf(1,1),!oro $ ftsol, !tsurf $ albsol, !albs $ pphi, !zma $ pphis, !phis $ cldfra, !cldfr $ rneb, !cldfr_st $ diafra, !cldfr_cv $ itop_con, !cldtop $ ibas_con, !cldbot $ cldliq, !cwat $ prfl, !flxrst $ pmflxr, !flxrcv $ psfl, !flxsst $ pmflxs, !flxscv $ pmfu, !flxupd $ flxmass_w, !flxmass_w $ t_seri, !tfld $ sh, !sh $ rh, !rh $ .false., !wrhstts $ hbuf, !hbuf $ obuf, !obuf $ iip1, !nx $ jjp1, !ny $ source, $ solsym) call VTe(VTinca) call VTb(VTphysiq) #endif END IF ! config_inca c====================================================================== c Calcul de l'effet de la convection c====================================================================== c print*,'Avant convection' do it=1,nqmax WRITE(itn,'(i2)') it c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn) enddo if (convection) then c print*,'Pas de temps dans phytrac : ',pdtphys DO it=1, nqmax IF ( config_inca/='none' .AND. conv_flg(it) == 0 ) CYCLE if (iflag_con.lt.2) then d_tr_cv=0. else if (iflag_con.eq.2) then c tiedke CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, . pplay, paprs, tr_seri(1,1,it), d_tr_cv(1,1,it)) else c KE call cvltr(pdtphys, da, phi, mp, paprs,pplay, tr_seri(1,1,it), . upwd,dnwd,d_tr_cv(1,1,it)) endif DO k = 1, nlev DO i = 1, klon tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k,it) ENDDO ENDDO IF (config_inca == 'none') THEN CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '//itn) ELSE CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = ' . //solsym(it)) END IF ENDDO endif ! convection c print*,'Apres convection' c do it=1,nqmax c WRITE(itn,'(i1)') it c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn) c enddo c====================================================================== c Calcul de l'effet des thermiques c====================================================================== do k=1,klev do i=1,klon zmasse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg enddo enddo c print*,'masse dans ph ',zmasse do it=1,nqmax do k=1,klev do i=1,klon d_tr_th(i,k,it)=0. tr_seri(i,k,it)=max(tr_seri(i,k,it),0.) tr_seri(i,k,it)=min(tr_seri(i,k,it),1.e10) enddo enddo enddo if (iflag_thermals.gt.0) then c print*,'calcul de leffet des thermiques' nsplit=10 DO it=1, nqmax c WRITE(itn,'(i1)') it c CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'conv it='//itn) c print*,'avant dqthermiquesretro' c call dump2d(iim,jjm-1,tr_seri(2,1,1),'TR_SERI ') do isplit=1,nsplit c Abderr 25 11 02 C Thermiques c print*,'Avant dans phytrac' call dqthermcell(klon,klev,pdtphys/nsplit . ,fm_therm,entr_therm,zmasse . ,tr_seri(1:klon,1:klev,it),d_tr,ztra_th) do k=1,klev do i=1,klon d_tr(i,k)=pdtphys*d_tr(i,k)/nsplit d_tr_th(i,k,it)=d_tr_th(i,k,it)+d_tr(i,k) tr_seri(i,k,it)=max(tr_seri(i,k,it)+d_tr(i,k),0.) enddo enddo enddo ! nsplit c print*,'apres thermiques' c call dump2d(iim,jjm-1,d_tr_th(1,1,1),'d_tr_th ') c do k=1,klev c print*,'d_tr_th(',k,')=',tr_seri(280,k,1) c enddo c WRITE(itn,'(i1)') it c CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'therm it='//itn) ENDDO ! it endif ! Thermiques c print*,'ATTENTION: sdans thermniques' c====================================================================== c Calcul de l'effet de la couche limite c====================================================================== c print *,'Avant couchelimite' c do it=1,nqmax c WRITE(itn,'(i1)') it c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn) c enddo if (couchelimite) then DO k = 1, nlev DO i = 1, klon delp(i,k) = paprs(i,k)-paprs(i,k+1) ENDDO ENDDO C maf modif pour tenir compte du cas rnpb + traceur DO it=1, nqmax IF ( config_inca/='none' .AND. pbl_flg(it) == 0 ) CYCLE c print *,'it',it,clsol(it) if (clsol(it)) then ! couche limite avec quantite dans le sol calculee CALL cltracrn(it, pdtphys, yu1, yv1, e cdragh, coefh,t_seri,ftsol,pctsrf, e tr_seri(1,1,it),trs(1,it), e paprs, pplay, delp, e masktr(1,it),fshtr(1,it),hsoltr(it), e tautr(it),vdeptr(it), e xlat, s d_tr_cl(1,1,it),d_trs) DO k = 1, nlev DO i = 1, klon tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it) ENDDO ENDDO c c Traceur ds sol c DO i = 1, klon trs(i,it) = trs(i,it) + d_trs(i) END DO C C maf provisoire suppression des prints C WRITE(itn,'(i1)') it C CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracrn it='//itn) else ! couche limite avec flux prescrit #ifndef INCA Cmaf provisoire source / traceur a creer DO i=1, klon source(i,it) = 0.0 ! pas de source, pour l'instant ENDDO C #endif CALL cltrac(pdtphys, coefh,t_seri, s tr_seri(1,1,it), source(:,it), e paprs, pplay, delp, s d_tr_cl(1,1,it)) DO k = 1, nlev DO i = 1, klon tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it) ENDDO ENDDO Cmaf WRITE(itn,'(i1)') it cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracn it='//itn) endif ENDDO c endif ! couche limite c print*,'Apres couchelimite' c do it=1,nqmax c WRITE(itn,'(i1)') it c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn) c enddo c====================================================================== c Calcul de l'effet du puits radioactif c====================================================================== C MAF il faudrait faire une modification pour passer dans radiornpb c si radio=true mais pour l'instant radiornpb propre au cas rnpb if(rnpb) then c print *, 'decroissance radiactive activee' call radiornpb (tr_seri,pdtphys,tautr,d_tr_dec) C DO it=1,nqmax if(radio(it)) then DO k = 1, nlev DO i = 1, klon tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it) ENDDO ENDDO WRITE(itn,'(i1)') it CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'puits rn it='//itn) endif ENDDO c endif ! rnpb decroissance radioactive C c====================================================================== c Calcul de l'effet de la precipitation c====================================================================== c print*,'LESSIVAGE =',lessivage IF (lessivage) THEN c print*,'avant lessivage' d_tr_lessi_nucl(:,:,:) = 0. d_tr_lessi_impa(:,:,:) = 0. flestottr(:,:,:) = 0. c c tendance des aerosols nuclees et impactes c DO it = 1, nqmax IF (aerosol(it)) THEN DO k = 1, nlev DO i = 1, klon d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) + s ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it) d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) + s ( 1 - frac_impa(i,k) )*tr_seri(i,k,it) ENDDO ENDDO ENDIF ENDDO c c Mises a jour des traceurs + calcul des flux de lessivage c Mise a jour due a l'impaction et a la nucleation c c call dump2d(iim,jjm-1,frac_impa(2:klon-1,10),'FRACIMPA') c call dump2d(iim,jjm-1,frac_nucl(2:klon-1,10),'FRACNUCL') c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3') DO it = 1, nqmax c print*,'IT=',it,aerosol(it) IF (aerosol(it)) THEN c print*,'IT=',it,' On lessive' DO k = 1, nlev DO i = 1, klon tr_seri(i,k,it)=tr_seri(i,k,it) s *frac_impa(i,k)*frac_nucl(i,k) ENDDO ENDDO ENDIF ENDDO c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3B') c c Flux lessivage total c DO it = 1, nqmax DO k = 1, nlev DO i = 1, klon flestottr(i,k,it) = flestottr(i,k,it) - s ( d_tr_lessi_nucl(i,k,it) + s d_tr_lessi_impa(i,k,it) ) * s ( paprs(i,k)-paprs(i,k+1) ) / s (RG * pdtphys) ENDDO ENDDO c Cmaf WRITE(itn,'(i1)') it Cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'tr(lessi) it='//itn) ENDDO c c print*,'apres lessivage' ENDIF Cc DO k = 1, nlev DO i = 1, klon fluxrn(i,k) = flestottr(i,k,1) fluxpb(i,k) = flestottr(i,k,2) rn(i,k) = tr_seri(i,k,1) pb(i,k) = tr_seri(i,k,2) pbnucl(i,k)=d_tr_lessi_nucl(i,k,2) pbimpa(i,k)=d_tr_lessi_impa(i,k,2) ENDDO ENDDO c============================================================= c Ecriture des sorties c============================================================= #ifdef CPP_IOIPSL #include "write_histrac.h" #endif c============================================================= if (lafin) then print*, 'c est la fin de la physique' call Gather(trs(:,1),trs_tmp) c$OMP MASTER if (is_mpi_root) then open (99,file='restarttrac', form='formatted') do i=1,klon_glo write(99,*) trs_tmp(i) enddo PRINT*, 'Ecriture du fichier restarttrac' close(99) endif c$OMP END MASTER else c print*, 'physique pas fini' endif RETURN END