! ! $Id: bilan_dyn_p.F 1550 2011-07-05 09:44:55Z acozic $ ! SUBROUTINE bilan_dyn_p (ntrac,dt_app,dt_cum, s ps,masse,pk,flux_u,flux_v,teta,phi,ucov,vcov,trac) c AFAIRE c Prevoir en champ nq+1 le diagnostique de l'energie c en faisant Qzon=Cv T + L * ... c vQ..A=Cp T + L * ... #ifdef CPP_IOIPSL USE IOIPSL #endif USE parallel USE mod_hallo use misc_mod use write_field IMPLICIT NONE #include "dimensions.h" #include "paramet.h" #include "comconst.h" #include "comvert.h" #include "comgeom2.h" #include "temps.h" #include "iniprint.h" c==================================================================== c c Sous-programme consacre à des diagnostics dynamiques de base c c c De facon generale, les moyennes des scalaires Q sont ponderees par c la masse. c c Les flux de masse sont eux simplement moyennes. c c==================================================================== c Arguments : c =========== integer ntrac real dt_app,dt_cum real ps(iip1,jjp1) real masse(iip1,jjp1,llm),pk(iip1,jjp1,llm) real flux_u(iip1,jjp1,llm) real flux_v(iip1,jjm,llm) real teta(iip1,jjp1,llm) real phi(iip1,jjp1,llm) real ucov(iip1,jjp1,llm) real vcov(iip1,jjm,llm) real trac(iip1,jjp1,llm,ntrac) c Local : c ======= integer icum,ncum logical first real zz,zqy,zfactv(jjm,llm) integer nQ parameter (nQ=7) cym character*6 nom(nQ) cym character*6 unites(nQ) character*6,save :: nom(nQ) character*6,save :: unites(nQ) character*10 file integer ifile parameter (ifile=4) integer itemp,igeop,iecin,iang,iu,iovap,iun integer i_sortie save first,icum,ncum save itemp,igeop,iecin,iang,iu,iovap,iun save i_sortie real time integer itau save time,itau data time,itau/0.,0/ data first/.true./ data itemp,igeop,iecin,iang,iu,iovap,iun/1,2,3,4,5,6,7/ data i_sortie/1/ real ww c variables dynamiques intermédiaires REAL vcont(iip1,jjm,llm),ucont(iip1,jjp1,llm) REAL ang(iip1,jjp1,llm),unat(iip1,jjp1,llm) REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm) REAL vorpot(iip1,jjm,llm) REAL w(iip1,jjp1,llm),ecin(iip1,jjp1,llm),convm(iip1,jjp1,llm) REAL bern(iip1,jjp1,llm) c champ contenant les scalaires advectés. real Q(iip1,jjp1,llm,nQ) c champs cumulés real ps_cum(iip1,jjp1) real masse_cum(iip1,jjp1,llm) real flux_u_cum(iip1,jjp1,llm) real flux_v_cum(iip1,jjm,llm) real Q_cum(iip1,jjp1,llm,nQ) real flux_uQ_cum(iip1,jjp1,llm,nQ) real flux_vQ_cum(iip1,jjm,llm,nQ) real flux_wQ_cum(iip1,jjp1,llm,nQ) real dQ(iip1,jjp1,llm,nQ) save ps_cum,masse_cum,flux_u_cum,flux_v_cum save Q_cum,flux_uQ_cum,flux_vQ_cum c champs de tansport en moyenne zonale integer ntr,itr parameter (ntr=5) cym character*10 znom(ntr,nQ) cym character*20 znoml(ntr,nQ) cym character*10 zunites(ntr,nQ) character*10,save :: znom(ntr,nQ) character*20,save :: znoml(ntr,nQ) character*10,save :: zunites(ntr,nQ) integer iave,itot,immc,itrs,istn data iave,itot,immc,itrs,istn/1,2,3,4,5/ character*3 ctrs(ntr) data ctrs/' ','TOT','MMC','TRS','STN'/ real zvQ(jjm,llm,ntr,nQ),zvQtmp(jjm,llm) real zavQ(jjm,ntr,nQ),psiQ(jjm,llm+1,nQ) real zmasse(jjm,llm),zamasse(jjm) real zv(jjm,llm),psi(jjm,llm+1) integer i,j,l,iQ c Initialisation du fichier contenant les moyennes zonales. c --------------------------------------------------------- character*10 infile integer fileid integer thoriid, zvertiid save fileid integer ndex3d(jjm*llm) C Variables locales C integer tau0 real zjulian character*3 str character*10 ctrac integer ii,jj integer zan, dayref C real rlong(jjm),rlatg(jjm) integer :: jjb,jje,jjn,ijb,ije type(Request) :: Req ! definition du domaine d'ecriture pour le rebuild INTEGER,DIMENSION(1) :: ddid INTEGER,DIMENSION(1) :: dsg INTEGER,DIMENSION(1) :: dsl INTEGER,DIMENSION(1) :: dpf INTEGER,DIMENSION(1) :: dpl INTEGER,DIMENSION(1) :: dhs INTEGER,DIMENSION(1) :: dhe INTEGER :: bilan_dyn_domain_id c===================================================================== c Initialisation c===================================================================== ndex3d=0 if (adjust) return time=time+dt_app itau=itau+1 if (first) then icum=0 c initialisation des fichiers first=.false. c ncum est la frequence de stokage en pas de temps ncum=dt_cum/dt_app if (abs(ncum*dt_app-dt_cum).gt.1.e-5*dt_app) then WRITE(lunout,*) . 'Pb : le pas de cumule doit etre multiple du pas' WRITE(lunout,*)'dt_app=',dt_app WRITE(lunout,*)'dt_cum=',dt_cum stop endif if (i_sortie.eq.1) then file='dynzon' if (mpi_rank==0) then call inigrads(ifile,1 s ,0.,180./pi,0.,0.,jjm,rlatv,-90.,90.,180./pi s ,llm,presnivs,1. s ,dt_cum,file,'dyn_zon ') endif endif nom(itemp)='T' nom(igeop)='gz' nom(iecin)='K' nom(iang)='ang' nom(iu)='u' nom(iovap)='ovap' nom(iun)='un' unites(itemp)='K' unites(igeop)='m2/s2' unites(iecin)='m2/s2' unites(iang)='ang' unites(iu)='m/s' unites(iovap)='kg/kg' unites(iun)='un' c Initialisation du fichier contenant les moyennes zonales. c --------------------------------------------------------- infile='dynzon' zan = annee_ref dayref = day_ref CALL ymds2ju(zan, 1, dayref, 0.0, zjulian) tau0 = itau_dyn rlong=0. rlatg=rlatv*180./pi jjb=jj_begin jje=jj_end jjn=jj_nb IF (pole_sud) THEN jjn=jj_nb-1 jje=jj_end-1 ENDIF ddid=(/ 2 /) dsg=(/ jjm /) dsl=(/ jjn /) dpf=(/ jjb /) dpl=(/ jje /) dhs=(/ 0 /) dhe=(/ 0 /) call flio_dom_set(mpi_size,mpi_rank,ddid,dsg,dsl,dpf,dpl,dhs,dhe, . 'box',bilan_dyn_domain_id) call histbeg(trim(infile), . 1, rlong(jjb:jje), jjn, rlatg(jjb:jje), . 1, 1, 1, jjn, . tau0, zjulian, dt_cum, thoriid, fileid, . bilan_dyn_domain_id) C C Appel a histvert pour la grille verticale C call histvert(fileid, 'presnivs', 'Niveaux sigma','mb', . llm, presnivs, zvertiid) C C Appels a histdef pour la definition des variables a sauvegarder do iQ=1,nQ do itr=1,ntr if(itr.eq.1) then znom(itr,iQ)=nom(iQ) znoml(itr,iQ)=nom(iQ) zunites(itr,iQ)=unites(iQ) else znom(itr,iQ)=ctrs(itr)//'v'//nom(iQ) znoml(itr,iQ)='transport : v * '//nom(iQ)//' '//ctrs(itr) zunites(itr,iQ)='m/s * '//unites(iQ) endif enddo enddo c Declarations des champs avec dimension verticale c print*,'1HISTDEF' do iQ=1,nQ do itr=1,ntr IF (prt_level > 5) . WRITE(lunout,*)'var ',itr,iQ . ,znom(itr,iQ),znoml(itr,iQ),zunites(itr,iQ) call histdef(fileid,znom(itr,iQ),znoml(itr,iQ), . zunites(itr,iQ),1,jjn,thoriid,llm,1,llm,zvertiid, . 32,'ave(X)',dt_cum,dt_cum) enddo c Declarations pour les fonctions de courant c print*,'2HISTDEF' call histdef(fileid,'psi'//nom(iQ) . ,'stream fn. '//znoml(itot,iQ), . zunites(itot,iQ),1,jjn,thoriid,llm,1,llm,zvertiid, . 32,'ave(X)',dt_cum,dt_cum) enddo c Declarations pour les champs de transport d'air c print*,'3HISTDEF' call histdef(fileid, 'masse', 'masse', . 'kg', 1, jjn, thoriid, llm, 1, llm, zvertiid, . 32, 'ave(X)', dt_cum, dt_cum) call histdef(fileid, 'v', 'v', . 'm/s', 1, jjn, thoriid, llm, 1, llm, zvertiid, . 32, 'ave(X)', dt_cum, dt_cum) c Declarations pour les fonctions de courant c print*,'4HISTDEF' call histdef(fileid,'psi','stream fn. MMC ','mega t/s', . 1,jjn,thoriid,llm,1,llm,zvertiid, . 32,'ave(X)',dt_cum,dt_cum) c Declaration des champs 1D de transport en latitude c print*,'5HISTDEF' do iQ=1,nQ do itr=2,ntr call histdef(fileid,'a'//znom(itr,iQ),znoml(itr,iQ), . zunites(itr,iQ),1,jjn,thoriid,1,1,1,-99, . 32,'ave(X)',dt_cum,dt_cum) enddo enddo c print*,'8HISTDEF' CALL histend(fileid) endif c===================================================================== c Calcul des champs dynamiques c ---------------------------- jjb=jj_begin jje=jj_end c énergie cinétique ucont(:,jjb:jje,:)=0 call Register_Hallo(ucov,ip1jmp1,llm,1,1,1,1,Req) call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Req) call SendRequest(Req) call WaitRequest(Req) CALL covcont_p(llm,ucov,vcov,ucont,vcont) CALL enercin_p(vcov,ucov,vcont,ucont,ecin) c moment cinétique do l=1,llm ang(:,jjb:jje,l)=ucov(:,jjb:jje,l)+constang(:,jjb:jje) unat(:,jjb:jje,l)=ucont(:,jjb:jje,l)*cu(:,jjb:jje) enddo Q(:,jjb:jje,:,itemp)=teta(:,jjb:jje,:)*pk(:,jjb:jje,:)/cpp Q(:,jjb:jje,:,igeop)=phi(:,jjb:jje,:) Q(:,jjb:jje,:,iecin)=ecin(:,jjb:jje,:) Q(:,jjb:jje,:,iang)=ang(:,jjb:jje,:) Q(:,jjb:jje,:,iu)=unat(:,jjb:jje,:) Q(:,jjb:jje,:,iovap)=trac(:,jjb:jje,:,1) Q(:,jjb:jje,:,iun)=1. c===================================================================== c Cumul c===================================================================== c if(icum.EQ.0) then jjb=jj_begin jje=jj_end ps_cum(:,jjb:jje)=0. masse_cum(:,jjb:jje,:)=0. flux_u_cum(:,jjb:jje,:)=0. Q_cum(:,jjb:jje,:,:)=0. flux_uQ_cum(:,jjb:jje,:,:)=0. if (pole_sud) jje=jj_end-1 flux_v_cum(:,jjb:jje,:)=0. flux_vQ_cum(:,jjb:jje,:,:)=0. endif IF (prt_level > 5) . WRITE(lunout,*)'dans bilan_dyn ',icum,'->',icum+1 icum=icum+1 c accumulation des flux de masse horizontaux jjb=jj_begin jje=jj_end ps_cum(:,jjb:jje)=ps_cum(:,jjb:jje)+ps(:,jjb:jje) masse_cum(:,jjb:jje,:)=masse_cum(:,jjb:jje,:)+masse(:,jjb:jje,:) flux_u_cum(:,jjb:jje,:)=flux_u_cum(:,jjb:jje,:) . +flux_u(:,jjb:jje,:) if (pole_sud) jje=jj_end-1 flux_v_cum(:,jjb:jje,:)=flux_v_cum(:,jjb:jje,:) . +flux_v(:,jjb:jje,:) jjb=jj_begin jje=jj_end do iQ=1,nQ Q_cum(:,jjb:jje,:,iQ)=Q_cum(:,jjb:jje,:,iQ) . +Q(:,jjb:jje,:,iQ)*masse(:,jjb:jje,:) enddo c===================================================================== c FLUX ET TENDANCES c===================================================================== c Flux longitudinal c ----------------- do iQ=1,nQ do l=1,llm do j=jjb,jje do i=1,iim flux_uQ_cum(i,j,l,iQ)=flux_uQ_cum(i,j,l,iQ) s +flux_u(i,j,l)*0.5*(Q(i,j,l,iQ)+Q(i+1,j,l,iQ)) enddo flux_uQ_cum(iip1,j,l,iQ)=flux_uQ_cum(1,j,l,iQ) enddo enddo enddo c flux méridien c ------------- do iQ=1,nQ call Register_Hallo(Q(1,1,1,iQ),ip1jmp1,llm,0,1,1,0,Req) enddo call SendRequest(Req) call WaitRequest(Req) jjb=jj_begin jje=jj_end if (pole_sud) jje=jj_end-1 do iQ=1,nQ do l=1,llm do j=jjb,jje do i=1,iip1 flux_vQ_cum(i,j,l,iQ)=flux_vQ_cum(i,j,l,iQ) s +flux_v(i,j,l)*0.5*(Q(i,j,l,iQ)+Q(i,j+1,l,iQ)) enddo enddo enddo enddo c tendances c --------- c convergence horizontale call Register_Hallo(flux_uQ_cum,ip1jmp1,llm,2,2,2,2,Req) call Register_Hallo(flux_vQ_cum,ip1jm,llm,2,2,2,2,Req) call SendRequest(Req) call WaitRequest(Req) call convflu_p(flux_uQ_cum,flux_vQ_cum,llm*nQ,dQ) c calcul de la vitesse verticale call Register_Hallo(flux_u_cum,ip1jmp1,llm,2,2,2,2,Req) call Register_Hallo(flux_v_cum,ip1jm,llm,2,2,2,2,Req) call SendRequest(Req) call WaitRequest(Req) call convmas_p(flux_u_cum,flux_v_cum,convm) CALL vitvert_p(convm,w) jjb=jj_begin jje=jj_end do iQ=1,nQ do l=1,llm-1 do j=jjb,jje do i=1,iip1 ww=-0.5*w(i,j,l+1)*(Q(i,j,l,iQ)+Q(i,j,l+1,iQ)) dQ(i,j,l ,iQ)=dQ(i,j,l ,iQ)-ww dQ(i,j,l+1,iQ)=dQ(i,j,l+1,iQ)+ww enddo enddo enddo enddo IF (prt_level > 5) . WRITE(lunout,*)'Apres les calculs fait a chaque pas' c===================================================================== c PAS DE TEMPS D'ECRITURE c===================================================================== if (icum.eq.ncum) then c===================================================================== IF (prt_level > 5) . WRITE(lunout,*)'Pas d ecriture' c Normalisation do iQ=1,nQ Q_cum(:,jjb:jje,:,iQ)=Q_cum(:,jjb:jje,:,iQ) . /masse_cum(:,jjb:jje,:) enddo zz=1./REAL(ncum) jjb=jj_begin jje=jj_end ps_cum(:,jjb:jje)=ps_cum(:,jjb:jje)*zz masse_cum(:,jjb:jje,:)=masse_cum(:,jjb:jje,:)*zz flux_u_cum(:,jjb:jje,:)=flux_u_cum(:,jjb:jje,:)*zz flux_uQ_cum(:,jjb:jje,:,:)=flux_uQ_cum(:,jjb:jje,:,:)*zz dQ(:,jjb:jje,:,:)=dQ(:,jjb:jje,:,:)*zz IF (pole_sud) jje=jj_end-1 flux_v_cum(:,jjb:jje,:)=flux_v_cum(:,jjb:jje,:)*zz flux_vQ_cum(:,jjb:jje,:,:)=flux_vQ_cum(:,jjb:jje,:,:)*zz jjb=jj_begin jje=jj_end c A retravailler eventuellement c division de dQ par la masse pour revenir aux bonnes grandeurs do iQ=1,nQ dQ(:,jjb:jje,:,iQ)=dQ(:,jjb:jje,:,iQ)/masse_cum(:,jjb:jje,:) enddo c===================================================================== c Transport méridien c===================================================================== c cumul zonal des masses des mailles c ---------------------------------- jjb=jj_begin jje=jj_end if (pole_sud) jje=jj_end-1 zv(jjb:jje,:)=0. zmasse(jjb:jje,:)=0. call Register_Hallo(masse_cum,ip1jmp1,llm,1,1,1,1,Req) do iQ=1,nQ call Register_Hallo(Q_cum(1,1,1,iQ),ip1jmp1,llm,0,1,1,0,Req) enddo call SendRequest(Req) call WaitRequest(Req) call massbar_p(masse_cum,massebx,masseby) jjb=jj_begin jje=jj_end if (pole_sud) jje=jj_end-1 do l=1,llm do j=jjb,jje do i=1,iim zmasse(j,l)=zmasse(j,l)+masseby(i,j,l) zv(j,l)=zv(j,l)+flux_v_cum(i,j,l) enddo zfactv(j,l)=cv(1,j)/zmasse(j,l) enddo enddo c print*,'3OK' c -------------------------------------------------------------- c calcul de la moyenne zonale du transport : c ------------------------------------------ c c -- c TOT : la circulation totale [ vq ] c c - - c MMC : mean meridional circulation [ v ] [ q ] c c ---- -- - - c TRS : transitoires [ v'q'] = [ vq ] - [ v q ] c c - * - * - - - - c STT : stationaires [ v q ] = [ v q ] - [ v ] [ q ] c c - - c on utilise aussi l'intermediaire TMP : [ v q ] c c la variable zfactv transforme un transport meridien cumule c en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte c c -------------------------------------------------------------- c ---------------------------------------- c Transport dans le plan latitude-altitude c ---------------------------------------- jjb=jj_begin jje=jj_end if (pole_sud) jje=jj_end-1 zvQ=0. psiQ=0. do iQ=1,nQ zvQtmp=0. do l=1,llm do j=jjb,jje c print*,'j,l,iQ=',j,l,iQ c Calcul des moyennes zonales du transort total et de zvQtmp do i=1,iim zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ) s +flux_vQ_cum(i,j,l,iQ) zqy= 0.5*(Q_cum(i,j,l,iQ)*masse_cum(i,j,l)+ s Q_cum(i,j+1,l,iQ)*masse_cum(i,j+1,l)) zvQtmp(j,l)=zvQtmp(j,l)+flux_v_cum(i,j,l)*zqy s /(0.5*(masse_cum(i,j,l)+masse_cum(i,j+1,l))) zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)+zqy enddo c print*,'aOK' c Decomposition zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)/zmasse(j,l) zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)*zfactv(j,l) zvQtmp(j,l)=zvQtmp(j,l)*zfactv(j,l) zvQ(j,l,immc,iQ)=zv(j,l)*zvQ(j,l,iave,iQ)*zfactv(j,l) zvQ(j,l,itrs,iQ)=zvQ(j,l,itot,iQ)-zvQtmp(j,l) zvQ(j,l,istn,iQ)=zvQtmp(j,l)-zvQ(j,l,immc,iQ) enddo enddo c fonction de courant meridienne pour la quantite Q do l=llm,1,-1 do j=jjb,jje psiQ(j,l,iQ)=psiQ(j,l+1,iQ)+zvQ(j,l,itot,iQ) enddo enddo enddo c fonction de courant pour la circulation meridienne moyenne psi(jjb:jje,:)=0. do l=llm,1,-1 do j=jjb,jje psi(j,l)=psi(j,l+1)+zv(j,l) zv(j,l)=zv(j,l)*zfactv(j,l) enddo enddo c print*,'4OK' c sorties proprement dites if (i_sortie.eq.1) then jjb=jj_begin jje=jj_end jjn=jj_nb if (pole_sud) jje=jj_end-1 if (pole_sud) jjn=jj_nb-1 do iQ=1,nQ do itr=1,ntr call histwrite(fileid,znom(itr,iQ),itau, s zvQ(jjb:jje,:,itr,iQ) s ,jjn*llm,ndex3d) enddo call histwrite(fileid,'psi'//nom(iQ), s itau,psiQ(jjb:jje,1:llm,iQ) s ,jjn*llm,ndex3d) enddo call histwrite(fileid,'masse',itau,zmasse(jjb:jje,1:llm) s ,jjn*llm,ndex3d) call histwrite(fileid,'v',itau,zv(jjb:jje,1:llm) s ,jjn*llm,ndex3d) psi(jjb:jje,:)=psi(jjb:jje,:)*1.e-9 call histwrite(fileid,'psi',itau,psi(jjb:jje,1:llm), s jjn*llm,ndex3d) endif c ----------------- c Moyenne verticale c ----------------- zamasse(jjb:jje)=0. do l=1,llm zamasse(jjb:jje)=zamasse(jjb:jje)+zmasse(jjb:jje,l) enddo zavQ(jjb:jje,:,:)=0. do iQ=1,nQ do itr=2,ntr do l=1,llm zavQ(jjb:jje,itr,iQ)=zavQ(jjb:jje,itr,iQ) s +zvQ(jjb:jje,l,itr,iQ) s *zmasse(jjb:jje,l) enddo zavQ(jjb:jje,itr,iQ)=zavQ(jjb:jje,itr,iQ)/zamasse(jjb:jje) call histwrite(fileid,'a'//znom(itr,iQ),itau, s zavQ(jjb:jje,itr,iQ),jjn*llm,ndex3d) enddo enddo c on doit pouvoir tracer systematiquement la fonction de courant. c===================================================================== c///////////////////////////////////////////////////////////////////// icum=0 !/////////////////////////////////////// endif ! icum.eq.ncum !/////////////////////////////////////// c///////////////////////////////////////////////////////////////////// c===================================================================== return end