c c $Header$ c SUBROUTINE phytrac (rnpb, I debutphy,lafin, I nqmax, I nlon,nlev,pdtphys, I t_seri,paprs,pplay, I pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, I coefh,yu1,yv1,ftsol,pctsrf,xlat, I frac_impa, frac_nucl, I xlon,presnivs,paire,pphis, O tr_seri) USE ioipsl USE histcom 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 "dimphy.h" #include "indicesol.h" #include "temps.h" #include "control.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 real pdtphys ! pas d'integration pour la physique (seconde) real t_seri(nlon,nlev) ! temperature real tr_seri(nlon,nlev,nbtr) ! traceur real paprs(nlon,nlev+1) ! pression pour chaque inter-couche (en Pa) real pplay(nlon,nlev) ! pression pour le mileu de chaque couche (en Pa) real presnivs(klev) ! pressions approximat. des milieux couches ( en PA) real paire(klon) real pphis(klon) logical debutphy ! le flag de l'initialisation de la physique logical lafin ! le flag de la fin de la physique integer ll c cAA Rem : nbtr : nombre de vrais traceurs est defini dans dimphy.h 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 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 c Couche limite: c -------------- c 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) ! a voir lorsque le flux est prescrit REAL topuits(klev,nbtr) ! a voir cAA cAA Pour la source de radon et son reservoir de sol cAA ................................................ REAL trs(klon,nbtr) ! Conc. radon ds le sol SAVE trs REAL masktr(klon,nbtr) ! Masque reservoir de sol traceur c Masque de l'echange avec la surface c (1 = reservoir) ou (possible => 1 ) SAVE masktr REAL fshtr(klon,nbtr) ! Flux surfacique dans le reservoir de sol SAVE fshtr REAL hsoltr(nbtr) ! Epaisseur equivalente du reservoir de sol SAVE hsoltr REAL tautr(nbtr) ! Constante de decroissance radioactive SAVE tautr REAL vdeptr(nbtr) ! Vitesse de depot sec dans la couche Brownienne SAVE vdeptr REAL scavtr(nbtr) ! Coefficient de lessivage SAVE scavtr cAA CHARACTER*1 itn C maf ioipsl CHARACTER*2 str2 INTEGER nhori, nvert REAL zsto, zout, zjulian INTEGER nid_tra SAVE nid_tra c REAL x(klon,klev,nbtr+2) ! traceurs INTEGER ndex(1) 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) INTEGER itra SAVE itra ! compteur pour la physique C C Variables liees a l'ecriture de la bande histoire : phytrac.nc c INTEGER ecrit_tra SAVE ecrit_tra logical ok_sync parameter (ok_sync = .true.) C C nature du traceur c logical aerosol(nbtr) ! Nature du traceur c ! aerosol(it) = true => aerosol c ! aerosol(it) = false => gaz c ! nat_trac(it) = 1. aerosolc logical clsol(nbtr) ! clsol(it) = true => CL sol calculee logical radio(nbtr) ! radio(it)=true => decroisssance radioactive save aerosol,clsol,radio C c====================================================================== c c Declaration des procedures appelees c c--modif convection tiedtke INTEGER i, k, it,itap save itap 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) ! tendance de traceurs couche limite REAL d_tr_cv(klon,klev) ! tendance de traceurs convection 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 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,sorties, s inirnpb data first,couchelimite,convection,lessivage,sorties c$$$ OM Test KE s /.true.,.true.,.true.,.true.,.true./ s /.true.,.true.,.false.,.true.,.true./ c c====================================================================== c ecrit_tra = NINT(86400./pdtphys *1.0) ! tous les jours modname='phytrac' c print*,'DANS PHYTRAC debutphy=',debutphy if (debutphy) then print*,'dans phytrac ',pdtphys,ecritphy,ecrit_tra ecrit_tra = NINT(86400./pdtphys/2.) ! tous les 12H c ecrit_tra = NINT(86400./pdtphys) ! tous les 24H if(nbtr.lt.nqmax) then print*,'NQMAX=',nqmax 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 itra=0 itap=0 C CALL ymds2ju(annee_ref, 1, 1, 0.0, zjulian) zjulian = zjulian + day_ref c CALL gr_fi_ecrit(1,klon,iim,jjm+1,xlon,zx_lon) DO i = 1, iim zx_lon(i,1) = xlon(i+1) zx_lon(i,jjm+1) = xlon(i+1) ENDDO c DO ll=1,klev c znivsig(ll)=float(ll) c ENDDO CALL gr_fi_ecrit(1,klon,iim,jjm+1,xlat,zx_lat) CALL histbeg("histrac", iim,zx_lon(:,1), jjm+1,zx_lat(1,:), . 1,iim,1,jjm+1, itau_phy, zjulian, pdtphys, . nhori, nid_tra) CALL histvert(nid_tra, "presnivs", "Vertical levels", "mb", . klev, presnivs, nvert) zsto = pdtphys zout = pdtphys * FLOAT(ecrit_tra) c CALL histdef(nid_tra, "phis", "Surface geop. height", "-", . iim,jjm+1,nhori, 1,1,1, -99, 32, . "once", zsto,zout) c CALL histdef(nid_tra, "aire", "Grid area", "-", . iim,jjm+1,nhori, 1,1,1, -99, 32, . "once", zsto,zout) goto 666 CALL histdef(nid_tra, "pyu1", "Vent niv 1", "-", . iim,jjm+1,nhori, 1,1,1, -99, 32, . "inst(X)", zsto,zout) CALL histdef(nid_tra, "pyv1", "Vent niv 1", "-", . iim,jjm+1,nhori, 1,1,1, -99, 32, . "inst(X)", zsto,zout) CALL histdef(nid_tra, "psrf1", "nature sol", "-", . iim,jjm+1,nhori, 1,1,1, -99, 32, . "inst(X)", zsto,zout) CALL histdef(nid_tra, "psrf2", "nature sol", "-", . iim,jjm+1,nhori, 1,1,1, -99, 32, . "inst(X)", zsto,zout) CALL histdef(nid_tra, "psrf3", "nature sol", "-", . iim,jjm+1,nhori, 1,1,1, -99, 32, . "inst(X)", zsto,zout) CALL histdef(nid_tra, "psrf4", "nature sol", "-", . iim,jjm+1,nhori, 1,1,1, -99, 32, . "inst(X)", zsto,zout) CALL histdef(nid_tra, "ftsol1", "temper sol", "-", . iim,jjm+1,nhori, 1,1,1, -99, 32, . "inst(X)", zsto,zout) CALL histdef(nid_tra, "ftsol2", "temper sol", "-", . iim,jjm+1,nhori, 1,1,1, -99, 32, . "inst(X)", zsto,zout) CALL histdef(nid_tra, "ftsol3", "temper sol", "-", . iim,jjm+1,nhori, 1,1,1, -99, 32, . "inst", zsto,zout) CALL histdef(nid_tra, "ftsol4", "temper sol", "-", . iim,jjm+1,nhori, 1,1,1, -99, 32, . "inst(X)", zsto,zout) CALL histdef(nid_tra, "pplay", "flux u mont","-", . iim,jjm+1,nhori, klev,1,klev,nvert, 32, . "inst(X)", zsto,zout) CALL histdef(nid_tra, "t", "flux u mont","-", . iim,jjm+1,nhori, klev,1,klev,nvert, 32, . "inst(X)", zsto,zout) CALL histdef(nid_tra, "mfu", "flux u mont","-", . iim,jjm+1,nhori, klev,1,klev,nvert, 32, . "ave(X)", zsto,zout) CALL histdef(nid_tra, "mfd", "flux u decen","-", . iim,jjm+1,nhori, klev,1,klev,nvert, 32, . "ave(X)", zsto,zout) CALL histdef(nid_tra, "en_u", "flux u mont","-", . iim,jjm+1,nhori, klev,1,klev,nvert, 32, . "ave(X)", zsto,zout) CALL histdef(nid_tra, "en_d", "flux u mont","-", . iim,jjm+1,nhori, klev,1,klev,nvert, 32, . "ave(X)", zsto,zout) CALL histdef(nid_tra, "de_u", "flux u mont","-", . iim,jjm+1,nhori, klev,1,klev,nvert, 32, . "ave(X)", zsto,zout) CALL histdef(nid_tra, "de_d", "flux u mont","-", . iim,jjm+1,nhori, klev,1,klev,nvert, 32, . "ave(X)", zsto,zout) CALL histdef(nid_tra, "coefh", "turbulent coef","-", . iim,jjm+1,nhori, klev,1,klev,nvert, 32, . "ave(X)", zsto,zout) 666 continue c DO it=1,nqmax IF (it.LE.99) THEN WRITE(str2,'(i2.2)') it C champ 3D CALL histdef(nid_tra, "tr"//str2, "Tracer No."//str2, "U/kga", . iim,jjm+1,nhori, klev,1,klev,nvert, 32, . "ave(X)", zsto,zout) IF (lessivage) THEN CALL histdef(nid_tra, "fl"//str2, "Tracer Flux No."//str2, . "U/m2/s",iim,jjm+1,nhori, klev,1,klev,nvert, 32, . "ave(X)", zsto,zout) ENDIF ELSE PRINT*, "Trop de traceurs" CALL abort ENDIF ENDDO CALL histend(nid_tra) ndex(1) = 0 c i = NINT(zout/zsto) CALL gr_fi_ecrit(1,klon,iim,jjm+1,pphis,zx_tmp_2d) CALL histwrite(nid_tra,"phis",i,zx_tmp_2d,iim*(jjm+1),ndex) C i = NINT(zout/zsto) CALL gr_fi_ecrit(1,klon,iim,jjm+1,paire,zx_tmp_2d) CALL histwrite(nid_tra,"aire",i,zx_tmp_2d,iim*(jjm+1),ndex) 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 do it=1,nqmax do i=1,klon trs(i,it) = 0. enddo END DO open (99,file='starttrac',status='old', . err=999,form='formatted') read(99,*) (trs(i,1),i=1,klon) 999 close(99) print*, 'apres starttrac' c Initialisation de la fraction d'aerosols lessivee c DO it=1,nqmax DO k = 1, nlev DO i = 1, klon d_tr_lessi_impa(i,k,it) = 0.0 d_tr_lessi_nucl(i,k,it) = 0.0 ENDDO ENDDO ENDDO 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(nqmax.gt.2) aerosol(3)=.true. c abder goto 777 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 ndex(1)=0 itap=itap+1 CALL gr_fi_ecrit(1,klon,iim,jjm+1,yu1,zx_tmp_2d) CALL histwrite(nid_tra,"pyu1",itap,zx_tmp_2d, s iim*(jjm+1),ndex) CALL gr_fi_ecrit(1,klon,iim,jjm+1,yv1,zx_tmp_2d) CALL histwrite(nid_tra,"pyv1",itap,zx_tmp_2d, s iim*(jjm+1),ndex) CALL gr_fi_ecrit(1,klon,iim,jjm+1,pftsol1,zx_tmp_2d) CALL histwrite(nid_tra,"ftsol1",itap,zx_tmp_2d, s iim*(jjm+1),ndex) CALL gr_fi_ecrit(1,klon,iim,jjm+1,pftsol2,zx_tmp_2d) CALL histwrite(nid_tra,"ftsol2",itap,zx_tmp_2d, s iim*(jjm+1),ndex) CALL gr_fi_ecrit(1,klon,iim,jjm+1,pftsol3,zx_tmp_2d) CALL histwrite(nid_tra,"ftsol3",itap,zx_tmp_2d, s iim*(jjm+1),ndex) CALL gr_fi_ecrit(1,klon,iim,jjm+1,pftsol4,zx_tmp_2d) CALL histwrite(nid_tra,"ftsol4",itap,zx_tmp_2d, s iim*(jjm+1),ndex) CALL gr_fi_ecrit(1,klon,iim,jjm+1,ppsrf1,zx_tmp_2d) CALL histwrite(nid_tra,"psrf1",itap,zx_tmp_2d, s iim*(jjm+1),ndex) CALL gr_fi_ecrit(1,klon,iim,jjm+1,ppsrf2,zx_tmp_2d) CALL histwrite(nid_tra,"psrf2",itap,zx_tmp_2d, s iim*(jjm+1),ndex) CALL gr_fi_ecrit(1,klon,iim,jjm+1,ppsrf3,zx_tmp_2d) CALL histwrite(nid_tra,"psrf3",itap,zx_tmp_2d, s iim*(jjm+1),ndex) CALL gr_fi_ecrit(1,klon,iim,jjm+1,ppsrf4,zx_tmp_2d) CALL histwrite(nid_tra,"psrf4",itap,zx_tmp_2d, s iim*(jjm+1),ndex) 777 continue c====================================================================== c Calcul de l'effet de la convection c====================================================================== print*,'Avant convection' do it=1,nqmax WRITE(itn,'(i1)') it c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn) enddo if (convection) then print*,'Pas de temps dans phytrac : ',pdtphys DO it=1, nqmax CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, . pplay, paprs, tr_seri(1,1,it), d_tr_cv) DO k = 1, nlev DO i = 1, klon tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k) ENDDO ENDDO c WRITE(itn,'(i1)') it c CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it='//itn) ENDDO c print*,'apres nflxtr' 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 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 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 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,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) 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 Cmaf provisoire source / traceur a creer DO i=1, klon source(i) = 0.0 ! pas de source, pour l'instant ENDDO C CALL cltrac(pdtphys, coefh,t_seri, s tr_seri(1,1,it), source, e paprs, pplay, delp, s d_tr ) DO k = 1, nlev DO i = 1, klon tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr(i,k) ENDDO ENDDO Cmaf provisoire suppression des prints 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 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====================================================================== print*,'LESSIVAGE =',lessivage IF (lessivage) THEN c print*,'avant lessivage' DO it = 1, nqmax DO k = 1, nlev DO i = 1, klon d_tr_lessi_nucl(i,k,it) = 0.0 d_tr_lessi_impa(i,k,it) = 0.0 flestottr(i,k,it) = 0.0 ENDDO ENDDO ENDDO 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 suppression provisoire 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 itra=itra+1 ndex(1) = 0 DO it=1,nqmax IF (it.LE.99) THEN WRITE(str2,'(i2.2)') it C champs 3D CALL gr_fi_ecrit(klev,klon,iim,jjm+1,tr_seri(1,1,it),zx_tmp_3d) CALL histwrite(nid_tra,"tr"//str2,itra,zx_tmp_3d, . iim*(jjm+1)*klev,ndex) c IF (lessivage) THEN c CALL gr_fi_ecrit(klev,klon,iim,jjm+1,flestottr(1,1,it),zx_tmp_3d) c CALL histwrite(nid_tra,"fl"//str2,itra,zx_tmp_3d, c . iim*(jjm+1)*klev,ndex) c ENDIF ELSE PRINT*, "Trop de traceurs" CALL abort ENDIF ENDDO goto 888 CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pplay,zx_tmp_3d) CALL histwrite(nid_tra,"pplay",itra,zx_tmp_3d, . iim*(jjm+1)*klev,ndex) CALL gr_fi_ecrit(klev,klon,iim,jjm+1,t_seri,zx_tmp_3d) CALL histwrite(nid_tra,"t",itra,zx_tmp_3d, . iim*(jjm+1)*klev,ndex) CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pmfu,zx_tmp_3d) CALL histwrite(nid_tra,"mfu",itra,zx_tmp_3d, . iim*(jjm+1)*klev,ndex) CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pmfd,zx_tmp_3d) CALL histwrite(nid_tra,"mfd",itra,zx_tmp_3d, . iim*(jjm+1)*klev,ndex) CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pen_u,zx_tmp_3d) CALL histwrite(nid_tra,"en_u",itra,zx_tmp_3d, . iim*(jjm+1)*klev,ndex) CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pen_d,zx_tmp_3d) CALL histwrite(nid_tra,"en_d",itra,zx_tmp_3d, . iim*(jjm+1)*klev,ndex) CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pde_d,zx_tmp_3d) CALL histwrite(nid_tra,"de_d",itra,zx_tmp_3d, . iim*(jjm+1)*klev,ndex) CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pde_u,zx_tmp_3d) CALL histwrite(nid_tra,"de_u",itra,zx_tmp_3d, . iim*(jjm+1)*klev,ndex) CALL gr_fi_ecrit(klev,klon,iim,jjm+1,coefh,zx_tmp_3d) CALL histwrite(nid_tra,"coefh",itra,zx_tmp_3d, . iim*(jjm+1)*klev,ndex) 888 continue c print*,'Sortie phytrac' c do it=1,nqmax c WRITE(itn,'(i1)') it c call diagtracphy(tr_seri(:,:,it),paprs,'Fin Phys '//itn) c enddo if (lafin) then print*, 'c est la fin de la physique' open (99,file='restarttrac', form='formatted') do i=1,klon write(99,*) trs(i,1) enddo PRINT*, 'Ecriture du fichier restarttrac' close(99) else print*, 'physique pas fini' endif RETURN END