Index: LMDZ5/trunk/libf/dyn3dmem/abort_gcm.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/abort_gcm.F (revision 1672) +++ LMDZ5/trunk/libf/dyn3dmem/abort_gcm.F (revision 1673) @@ -1,4 +1,4 @@ ! -! $Id: abort_gcm.F 1425 2010-09-02 13:45:23Z lguez $ +! $Id$ ! c @@ -45,5 +45,4 @@ if (ierr .eq. 0) then write(lunout,*) 'Everything is cool' - stop else write(lunout,*) 'Houston, we have a problem ', ierr Index: LMDZ5/trunk/libf/dyn3dmem/academic.h =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/academic.h (revision 1672) +++ LMDZ5/trunk/libf/dyn3dmem/academic.h (revision 1673) @@ -1,5 +1,9 @@ ! -! $Header$ +! $Id$ ! - real tetarappel(ip1jmp1,llm),taurappel - common/academic/tetarappel,taurappel + common/academic/tetarappel,knewt_t,kfrict,knewt_g,clat4 + real :: tetarappel(ip1jmp1,llm) + real :: knewt_t(llm) + real :: kfrict(llm) + real :: knewt_g + real :: clat4(ip1jmp1) Index: LMDZ5/trunk/libf/dyn3dmem/addfi_loc.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/addfi_loc.F (revision 1672) +++ LMDZ5/trunk/libf/dyn3dmem/addfi_loc.F (revision 1673) @@ -1,4 +1,4 @@ ! -! $Header$ +! $Id$ ! SUBROUTINE addfi_loc(pdt, leapf, forward, @@ -7,4 +7,5 @@ USE parallel USE infotrac, ONLY : nqtot + USE control_mod, ONLY : planet_type IMPLICIT NONE c @@ -156,4 +157,6 @@ c$OMP END MASTER + if (planet_type=="earth") then + ! earth case, special treatment for first 2 tracers (water) DO iq = 1, 2 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) @@ -177,4 +180,17 @@ c$OMP END DO NOWAIT ENDDO + else + ! general case, treat all tracers equally) + DO iq = 1, nqtot +c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) + DO k = 1,llm + DO j = ijb,ije + pq(j,k,iq)= pq(j,k,iq) + pdqfi(j,k,iq) * pdt + pq(j,k,iq)= AMAX1( pq(j,k,iq), qtestt ) + ENDDO + ENDDO +c$OMP END DO NOWAIT + ENDDO + endif ! of if (planet_type=="earth") c$OMP MASTER Index: LMDZ5/trunk/libf/dyn3dmem/advtrac_loc.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/advtrac_loc.F (revision 1672) +++ LMDZ5/trunk/libf/dyn3dmem/advtrac_loc.F (revision 1673) @@ -1,4 +1,4 @@ ! -! $Id: advtrac_p.F 1299 2010-01-20 14:27:21Z fairhead $ +! $Id$ ! c @@ -342,4 +342,6 @@ c$OMP BARRIER + if (planet_type=="earth") then + ijb=ij_begin ije=ij_end @@ -355,4 +357,5 @@ CALL qminimum_loc( q, 2, finmasse ) + endif ! of if (planet_type=="earth") RETURN Index: LMDZ5/trunk/libf/dyn3dmem/bands.F90 =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/bands.F90 (revision 1672) +++ LMDZ5/trunk/libf/dyn3dmem/bands.F90 (revision 1673) @@ -1,4 +1,4 @@ ! -! $Id: bands.F90 1279 2009-12-10 09:02:56Z fairhead $ +! $Id$ ! module Bands @@ -105,6 +105,6 @@ SUBROUTINE Set_Bands USE parallel -#ifdef CPP_EARTH -! Ehouarn: what follows is only related to // physics; for now only for Earth +#ifdef CPP_PHYS +! Ehouarn: what follows is only related to // physics USE mod_phys_lmdz_para, ONLY : jj_para_begin,jj_para_end #endif @@ -118,6 +118,5 @@ enddo -#ifdef CPP_EARTH -! Ehouarn: what follows is only related to // physics; for now only for Earth +#ifdef CPP_PHYS do i=0,MPI_Size-1 jj_Nb_physic(i)=jj_para_end(i)-jj_para_begin(i)+1 @@ -374,6 +373,6 @@ subroutine AdjustBands_physic use times -#ifdef CPP_EARTH -! Ehouarn: what follows is only related to // physics; for now only for Earth +#ifdef CPP_PHYS +! Ehouarn: what follows is only related to // physics USE mod_phys_lmdz_para, only : klon_mpi_para_nb #endif @@ -401,6 +400,5 @@ medium=medium/mpi_size NbTot=0 -#ifdef CPP_EARTH -! Ehouarn: what follows is only related to // physics; for now only for Earth +#ifdef CPP_PHYS do i=0,mpi_size-1 Inc(i)=nint(klon_mpi_para_nb(i)*(medium-value(i))/value(i)) Index: LMDZ5/trunk/libf/dyn3dmem/bilan_dyn_loc.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/bilan_dyn_loc.F (revision 1672) +++ LMDZ5/trunk/libf/dyn3dmem/bilan_dyn_loc.F (revision 1673) @@ -421,5 +421,4 @@ Q_cum(:,jjb:jje,:,l)=0. flux_uQ_cum(:,jjb:jje,l,:)=0. - flux_v_cum(:,jjb:jje,l)=0. if (pole_sud) jje=jj_end-1 flux_v_cum(:,jjb:jje,l)=0. Index: LMDZ5/trunk/libf/dyn3dmem/bilan_dyn_p.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/bilan_dyn_p.F (revision 1672) +++ (revision ) @@ -1,717 +1,0 @@ -! -! $Id: bilan_dyn_p.F 1403 2010-07-01 09:02:53Z fairhead $ -! - 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. - flux_v_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 Index: LMDZ5/trunk/libf/dyn3dmem/caladvtrac_loc.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/caladvtrac_loc.F (revision 1672) +++ LMDZ5/trunk/libf/dyn3dmem/caladvtrac_loc.F (revision 1673) @@ -8,6 +8,6 @@ * flxw, pk, iapptrac) USE parallel - USE infotrac - USE control_mod + USE infotrac, ONLY : nqtot + USE control_mod, ONLY : iapp_tracvl,planet_type USE caladvtrac_mod USE mod_hallo @@ -38,5 +38,5 @@ REAL :: masse(ijb_u:ije_u,llm) REAL :: p( ijb_u:ije_u,llmp1) - REAL :: q( ijb_u:ije_u,llm,nqtot),dq( ijb_u:ije_u,llm,2 ) + REAL :: q( ijb_u:ije_u,llm,nqtot),dq( ijb_u:ije_u,llm, nqtot ) REAL :: teta( ijb_u:ije_u,llm),pk( ijb_u:ije_u,llm) REAL :: flxw(ijb_u:ije_u,llm) Index: LMDZ5/trunk/libf/dyn3dmem/calfis_loc.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/calfis_loc.F (revision 1672) +++ LMDZ5/trunk/libf/dyn3dmem/calfis_loc.F (revision 1673) @@ -27,6 +27,6 @@ $ pdqfi, $ pdpsfi) -#ifdef CPP_EARTH -! Ehouarn: For now, calfis_p needs Earth physics +#ifdef CPP_PHYS +! If using physics c c Auteur : P. Le Van, F. Hourdin @@ -36,8 +36,9 @@ USE parallel, ONLY : omp_chunk, using_mpi,jjb_u,jje_u,jjb_v,jje_v USE mod_interface_dyn_phys + USE IOPHY +#endif USE Write_Field Use Write_field_p USE Times - USE IOPHY USE infotrac USE control_mod @@ -145,4 +146,6 @@ +#ifdef CPP_PHYS +! Ehouarn: for now calfis_p needs some informations from physics to compile c Local variables : c ----------------- @@ -220,5 +223,5 @@ PARAMETER(ntetaSTD=3) REAL rtetaSTD(ntetaSTD) - DATA rtetaSTD/350., 380., 405./ + DATA rtetaSTD/350., 380., 405./ ! Earth-specific values, beware !! REAL PVteta(klon,ntetaSTD) @@ -243,5 +246,7 @@ REAL,ALLOCATABLE,SAVE:: zdufi2(:,:),zdvfi2(:,:) integer :: k,kstart,kend - INTEGER :: offset + INTEGER :: offset + + LOGICAL tracerdyn c c----------------------------------------------------------------------- @@ -512,6 +517,7 @@ - IF (is_sequential) THEN -c + IF (is_sequential.and.(planet_type=="earth")) THEN +#ifdef CPP_PHYS +! PVtheta calls tetalevel, which is in the physics cIM calcul PV a teta=350, 380, 405K CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta, @@ -519,4 +525,5 @@ $ ntetaSTD,rtetaSTD,PVteta) c +#endif ENDIF @@ -662,7 +669,4 @@ c$OMP BARRIER - if (planet_type=="earth") then -#ifdef CPP_EARTH - !$OMP MASTER @@ -675,4 +679,5 @@ zdqfic_omp(:,:,:)=0. +#ifdef CPP_PHYS do isplit=1,nsplit_phys @@ -681,4 +686,5 @@ lafin_split=lafin.and.isplit==nsplit_phys + if (planet_type=="earth") then CALL physiq (klon, @@ -711,4 +717,33 @@ . PVteta) + else if ( planet_type=="generic" ) then + + CALL physiq (klon, !! ngrid + . llm, !! nlayer + . nqtot, !! nq + . tname, !! tracer names from dynamical core (given in infotrac) + . debut_split, !! firstcall + . lafin_split, !! lastcall + . float(day_ini), !! pday <-- day_ini (dans temps.h) + . jH_cur_split, !! ptime "fraction of day" + . zdt_split, !! ptimestep + . zplev_omp, !! pplev + . zplay_omp, !! pplay + . zphi_omp, !! pphi + . zufi_omp, !! pu + . zvfi_omp, !! pv + . ztfi_omp, !! pt + . zqfi_omp, !! pq + . flxwfi_omp, !! pw !! or 0. anyway this is for diagnostic. not used in physiq. + . zdufi_omp, !! pdu + . zdvfi_omp, !! pdv + . zdtfi_omp, !! pdt + . zdqfi_omp, !! pdq + . zdpsrf_omp, !! pdpsrf + . tracerdyn) !! tracerdyn <-- utilite ??? + + endif ! of if (planet_type=="earth") + + zufi_omp(:,:)=zufi_omp(:,:)+zdufi_omp(:,:)*zdt_split zvfi_omp(:,:)=zvfi_omp(:,:)+zdvfi_omp(:,:)*zdt_split @@ -723,4 +758,8 @@ enddo +#endif +! of #ifdef CPP_PHYS + + zdufi_omp(:,:)=zdufic_omp(:,:)/nsplit_phys zdvfi_omp(:,:)=zdvfic_omp(:,:)/nsplit_phys @@ -728,6 +767,4 @@ zdqfi_omp(:,:,:)=zdqfic_omp(:,:,:)/nsplit_phys -#endif - endif !of if (planet_type=="earth") c$OMP BARRIER @@ -1179,11 +1216,10 @@ firstcal = .FALSE. -#else - write(*,*) "calfis_p: for now can only work with parallel physics" - write(lunout,*) - & "calfis_p: for now can only work with parallel physics" - stop -#endif -! of #ifdef CPP_EARTH +#else + write(lunout,*) + & "calfis_p: for now can only work with parallel physics" + stop +#endif +! of #ifdef CPP_PHYS RETURN END Index: LMDZ5/trunk/libf/dyn3dmem/call_calfis_mod.F90 =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/call_calfis_mod.F90 (revision 1672) +++ LMDZ5/trunk/libf/dyn3dmem/call_calfis_mod.F90 (revision 1673) @@ -136,6 +136,12 @@ !$OMP END MASTER - jD_cur = jD_ref + day_ini - day_ref + int (itau * dtvr / daysec) - jH_cur = jH_ref + (itau * dtvr / daysec - int(itau * dtvr / daysec)) + jD_cur = jD_ref + day_ini - day_ref & + & + itau/day_step + jH_cur = jH_ref + start_time + & + & mod(itau,day_step)/float(day_step) + if (jH_cur > 1.0 ) then + jD_cur = jD_cur +1. + jH_cur = jH_cur -1. + endif ! Inbterface avec les routines de phylmd (phymars ... ) Index: LMDZ5/trunk/libf/dyn3dmem/call_dissip_mod.F90 =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/call_dissip_mod.F90 (revision 1672) +++ LMDZ5/trunk/libf/dyn3dmem/call_dissip_mod.F90 (revision 1673) @@ -240,4 +240,8 @@ !$OMP END DO NOWAIT + if (1 == 0) then +!!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics +!!! 2) should probably not be here anyway +!!! but are kept for those who would want to revert to previous behaviour !$OMP MASTER DO ij = 1,iim @@ -251,5 +255,6 @@ !$OMP END MASTER - ENDIF + ENDIF ! of if (1 == 0) + endif ! of of (pole_nord) IF (pole_sud) THEN @@ -269,4 +274,8 @@ !$OMP END DO NOWAIT + if (1 == 0) then +!!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics +!!! 2) should probably not be here anyway +!!! but are kept for those who would want to revert to previous behaviour !$OMP MASTER DO ij = 1,iim @@ -279,5 +288,6 @@ ENDDO !$OMP END MASTER - ENDIF + ENDIF ! of if (1 == 0) + endif ! of if (pole_sud) Index: LMDZ5/trunk/libf/dyn3dmem/ce0l.F90 =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/ce0l.F90 (revision 1672) +++ LMDZ5/trunk/libf/dyn3dmem/ce0l.F90 (revision 1673) @@ -1,4 +1,4 @@ ! -! $Id: ce0l.F90 1425 2010-09-02 13:45:23Z lguez $ +! $Id$ ! !------------------------------------------------------------------------------- @@ -19,7 +19,8 @@ USE dimphy USE comgeomphy - USE mod_phys_lmdz_para + USE mod_phys_lmdz_para, mpi_root_xx=>mpi_root USE mod_const_mpi USE infotrac + USE parallel, ONLY: finalize_parallel #ifdef CPP_IOIPSL @@ -30,4 +31,5 @@ IMPLICIT NONE #ifndef CPP_EARTH +#include "iniprint.h" WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics' #else @@ -41,12 +43,20 @@ #include "temps.h" #include "logic.h" +#ifdef CPP_MPI + include 'mpif.h' +#endif + INTEGER, PARAMETER :: longcles=20 + INTEGER :: ierr REAL, DIMENSION(longcles) :: clesphy0 REAL, DIMENSION(iip1,jjp1) :: masque CHARACTER(LEN=15) :: calnd + REAL, DIMENSION(iip1,jjp1) :: phis ! geopotentiel au sol !------------------------------------------------------------------------------- CALL conf_gcm( 99, .TRUE. , clesphy0 ) +#ifdef CPP_MPI CALL init_mpi +#endif CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/)) @@ -55,5 +65,5 @@ CALL abort_gcm('ce0l','In parallel mode, & & ce0l must be called only & - & for 1 process and 1 task') + & for 1 process and 1 task',1) ENDIF @@ -76,9 +86,12 @@ #endif - IF (config_inca /= 'none') THEN + IF (type_trac == 'inca') THEN #ifdef INCA - CALL init_const_lmdz(nbtr,anneeref,dayref,iphysiq,day_step,nday) - CALL init_inca_para(iim,jjm+1,klon,1,klon_mpi_para_nb,0) - WRITE(lunout,*)'nbtr =' , nbtr + CALL init_const_lmdz( & + nbtr,anneeref,dayref,& + iphysiq,day_step,nday,& + nbsrf, is_oce,is_sic,& + is_ter,is_lic) + #endif END IF @@ -90,5 +103,5 @@ WRITE(lunout,'(//)') WRITE(lunout,*) ' interbar = ',interbar - CALL etat0_netcdf(interbar,masque,ok_etat0) + CALL etat0_netcdf(interbar,masque,phis,ok_etat0) IF(ok_limit) THEN @@ -101,4 +114,20 @@ END IF + IF (grilles_gcm_netcdf) THEN + WRITE(lunout,'(//)') + WRITE(lunout,*) ' *************************** ' + WRITE(lunout,*) ' *** grilles_gcm_netcdf *** ' + WRITE(lunout,*) ' *************************** ' + WRITE(lunout,'(//)') + CALL grilles_gcm_netcdf_sub(masque,phis) + END IF + +#ifdef CPP_MPI +!$OMP MASTER + CALL MPI_FINALIZE(ierr) + IF (ierr /= 0) CALL abort_gcm('ce0l','Error in MPI_FINALIZE',1) +!$OMP END MASTER +#endif + #endif ! of #ifndef CPP_EARTH #else Index: LMDZ5/trunk/libf/dyn3dmem/comconst.h =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/comconst.h (revision 1672) +++ LMDZ5/trunk/libf/dyn3dmem/comconst.h (revision 1673) @@ -1,23 +1,40 @@ ! -! $Id: comconst.h 1279 2009-12-10 09:02:56Z fairhead $ +! $Id$ ! !----------------------------------------------------------------------- ! INCLUDE comconst.h - COMMON/comconst/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl, & - & dtvr,daysec, & + COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl, & + & iflag_top_bound + COMMON/comconstr/dtvr,daysec, & & pi,dtphys,dtdiss,rad,r,cpp,kappa,cotot,unsim,g,omeg & & ,dissip_factz,dissip_deltaz,dissip_zref & - & ,iflag_top_bound,tau_top_bound + & ,tau_top_bound, & + & daylen,year_day,molmass, ihf INTEGER im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl - REAL dtvr,daysec - REAL pi,dtphys,dtdiss,rad,r,cpp,kappa - REAL cotot,unsim,g,omeg + REAL dtvr ! dynamical time step (in s) + REAL daysec !length (in s) of a standard day + REAL pi ! something like 3.14159.... + REAL dtphys ! (s) time step for the physics + REAL dtdiss ! (s) time step for the dissipation + REAL rad ! (m) radius of the planet + REAL r ! Reduced Gas constant r=R/mu + ! with R=8.31.. J.K-1.mol-1, mu: mol mass of atmosphere (kg/mol) + REAL cpp ! Specific heat Cp (J.kg-1.K-1) + REAL kappa ! kappa=R/Cp + REAL cotot + REAL unsim ! = 1./iim + REAL g ! (m/s2) gravity + REAL omeg ! (rad/s) rotation rate of the planet REAL dissip_factz,dissip_deltaz,dissip_zref INTEGER iflag_top_bound REAL tau_top_bound + REAL daylen ! length of solar day, in 'standard' day length + REAL year_day ! Number of standard days in a year + REAL molmass ! (g/mol) molar mass of the atmosphere + REAL ihf ! (W/m2) Intrinsic heat flux (for giant planets) !----------------------------------------------------------------------- Index: LMDZ5/trunk/libf/dyn3dmem/comdissipn.h =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/comdissipn.h (revision 1672) +++ LMDZ5/trunk/libf/dyn3dmem/comdissipn.h (revision 1673) @@ -1,16 +1,20 @@ ! -! $Header$ +! $Id$ ! -c----------------------------------------------------------------------- -c INCLUDE comdissipn.h +! Attention : ce fichier include est compatible format fixe/format libre +! veillez à n'utiliser que des ! pour les commentaires +! et à bien positionner les & des lignes de continuation +! (les placer en colonne 6 et en colonne 73) +!----------------------------------------------------------------------- +! INCLUDE comdissipn.h REAL tetaudiv, tetaurot, tetah, cdivu, crot, cdivh -c - COMMON/comdissipn/ tetaudiv(llm),tetaurot(llm),tetah(llm) , - 1 cdivu, crot, cdivh +! + COMMON/comdissipn/ tetaudiv(llm),tetaurot(llm),tetah(llm) , & + & cdivu, crot, cdivh -c -c Les parametres de ce common proviennent des calculs effectues dans -c Inidissip . -c -c----------------------------------------------------------------------- +! +! Les parametres de ce common proviennent des calculs effectues dans +! Inidissip . +! +!----------------------------------------------------------------------- Index: LMDZ5/trunk/libf/dyn3dmem/comvert.h =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/comvert.h (revision 1672) +++ LMDZ5/trunk/libf/dyn3dmem/comvert.h (revision 1673) @@ -1,14 +1,27 @@ ! -! $Header$ +! $Id$ ! !----------------------------------------------------------------------- ! INCLUDE 'comvert.h' - COMMON/comvert/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm), & - & pa,preff,nivsigs(llm),nivsig(llm+1) + COMMON/comvertr/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm), & + & pa,preff,nivsigs(llm),nivsig(llm+1), & + & aps(llm),bps(llm),scaleheight,pseudoalt(llm) - common/comverti/disvert_type + common/comverti/disvert_type, pressure_exner - REAL ap,bp,presnivs,dpres,pa,preff,nivsigs,nivsig + real ap ! hybrid pressure contribution at interlayers + real bp ! hybrid sigma contribution at interlayer + real presnivs ! (reference) pressure at mid-layers + real dpres + real pa ! reference pressure (Pa) at which hybrid coordinates + ! become purely pressure + real preff ! reference surface pressure (Pa) + real nivsigs + real nivsig + real aps ! hybrid pressure contribution at mid-layers + real bps ! hybrid sigma contribution at mid-layers + real scaleheight ! atmospheric (reference) scale height (km) + real pseudoalt ! for planets integer disvert_type ! type of vertical discretization: @@ -17,3 +30,8 @@ ! 2: Planets (default for planet_type!=earth), ! using 'z2sig.def' (or 'esasig.def) file -!----------------------------------------------------------------------- + + logical pressure_exner +! compute pressure inside layers using Exner function, else use mean +! of pressure values at interfaces + + !----------------------------------------------------------------------- Index: LMDZ5/trunk/libf/dyn3dmem/conf_gcm.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/conf_gcm.F (revision 1672) +++ LMDZ5/trunk/libf/dyn3dmem/conf_gcm.F (revision 1673) @@ -1,4 +1,4 @@ ! -! $Id: conf_gcm.F 1403 2010-07-01 09:02:53Z fairhead $ +! $Id$ ! c @@ -18,4 +18,5 @@ use parallel, ONLY : omp_chunk USE control_mod + USE infotrac, ONLY : type_trac IMPLICIT NONE c----------------------------------------------------------------------- @@ -103,5 +104,7 @@ CALL getin('lunout', lunout) IF (lunout /= 5 .and. lunout /= 6) THEN - OPEN(lunout,FILE='lmdz.out') + OPEN(UNIT=lunout,FILE='lmdz.out_0000',ACTION='write', + & STATUS='unknown',FORM='formatted') + ENDIF @@ -166,4 +169,12 @@ CALL getin('nday',nday) +!Config Key = starttime +!Config Desc = Heure de depart de la simulation +!Config Def = 0 +!Config Help = Heure de depart de la simulation +!Config en jour + starttime = 0 + CALL getin('starttime',starttime) + !Config Key = day_step !Config Desc = nombre de pas par jour @@ -226,11 +237,12 @@ CALL getin('output_grads_dyn',output_grads_dyn) -!Config Key = idissip +!Config Key = dissip_period !Config Desc = periode de la dissipation -!Config Def = 10 +!Config Def = 0 !Config Help = periode de la dissipation -!Config (en pas) ... a completer ! - idissip = 10 - CALL getin('idissip',idissip) +!Config dissip_period=0 => la valeur sera calcule dans inidissip +!Config dissip_period>0 => on prend cette valeur + dissip_period = 0 + CALL getin('dissip_period',dissip_period) ccc .... P. Le Van , modif le 29/04/97 .pour la dissipation ... @@ -579,4 +591,23 @@ offline = .FALSE. CALL getin('offline',offline) + IF (offline .AND. adjust) THEN + WRITE(lunout,*) + & 'WARNING : option offline does not work with adjust=y :' + WRITE(lunout,*) 'the files defstoke.nc, fluxstoke.nc ', + & 'and fluxstokev.nc will not be created' + WRITE(lunout,*) + & 'only the file phystoke.nc will still be created ' + END IF + +!Config Key = type_trac +!Config Desc = Choix de couplage avec model de chimie INCA ou REPROBUS +!Config Def = lmdz +!Config Help = +!Config 'lmdz' = pas de couplage, pur LMDZ +!Config 'inca' = model de chime INCA +!Config 'repr' = model de chime REPROBUS + type_trac = 'lmdz' + CALL getin('type_trac',type_trac) + !Config Key = config_inca @@ -628,5 +659,5 @@ write(lunout,*)' periodav = ', periodav write(lunout,*)' output_grads_dyn = ', output_grads_dyn - write(lunout,*)' idissip = ', idissip + write(lunout,*)' dissip_period = ', dissip_period write(lunout,*)' lstardis = ', lstardis write(lunout,*)' nitergdiv = ', nitergdiv @@ -651,4 +682,5 @@ write(lunout,*)' tauyy = ', tauyy write(lunout,*)' offline = ', offline + write(lunout,*)' type_trac = ', type_trac write(lunout,*)' config_inca = ', config_inca write(lunout,*)' ok_dynzon = ', ok_dynzon @@ -769,4 +801,22 @@ offline = .FALSE. CALL getin('offline',offline) + IF (offline .AND. adjust) THEN + WRITE(lunout,*) + & 'WARNING : option offline does not work with adjust=y :' + WRITE(lunout,*) 'the files defstoke.nc, fluxstoke.nc ', + & 'and fluxstokev.nc will not be created' + WRITE(lunout,*) + & 'only the file phystoke.nc will still be created ' + END IF + +!Config Key = type_trac +!Config Desc = Choix de couplage avec model de chimie INCA ou REPROBUS +!Config Def = lmdz +!Config Help = +!Config 'lmdz' = pas de couplage, pur LMDZ +!Config 'inca' = model de chime INCA +!Config 'repr' = model de chime REPROBUS + type_trac = 'lmdz' + CALL getin('type_trac',type_trac) !Config Key = config_inca @@ -875,4 +925,10 @@ CALL getin('ok_etat0',ok_etat0) +!Config Key = grilles_gcm_netcdf +!Config Desc = creation de fichier grilles_gcm.nc dans create_etat0_limit +!Config Def = n + grilles_gcm_netcdf = .FALSE. + CALL getin('grilles_gcm_netcdf',grilles_gcm_netcdf) + write(lunout,*)' #########################################' write(lunout,*)' Configuration des parametres de cel0' @@ -889,5 +945,5 @@ write(lunout,*)' periodav = ', periodav write(lunout,*)' output_grads_dyn = ', output_grads_dyn - write(lunout,*)' idissip = ', idissip + write(lunout,*)' dissip_period = ', dissip_period write(lunout,*)' lstardis = ', lstardis write(lunout,*)' nitergdiv = ', nitergdiv @@ -912,4 +968,5 @@ write(lunout,*)' tauy = ', tauy write(lunout,*)' offline = ', offline + write(lunout,*)' type_trac = ', type_trac write(lunout,*)' config_inca = ', config_inca write(lunout,*)' ok_dynzon = ', ok_dynzon @@ -923,4 +980,5 @@ write(lunout,*)' ok_limit = ', ok_limit write(lunout,*)' ok_etat0 = ', ok_etat0 + write(lunout,*)' grilles_gcm_netcdf = ', grilles_gcm_netcdf c RETURN Index: LMDZ5/trunk/libf/dyn3dmem/conf_planete.F90 =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/conf_planete.F90 (revision 1673) +++ LMDZ5/trunk/libf/dyn3dmem/conf_planete.F90 (revision 1673) @@ -0,0 +1,70 @@ +! +! $Id$ +! +SUBROUTINE conf_planete +! +#ifdef CPP_IOIPSL +USE IOIPSL +#else +! if not using IOIPSL, we still need to use (a local version of) getin +USE ioipsl_getincom +#endif +IMPLICIT NONE +! +! +! Declarations : +! -------------- +#include "dimensions.h" +#include "comconst.h" +#include "comvert.h" +! +! local: +! ------ + +! --------------------------------------------- +! Initialisations de constantes de la dynamique +! --------------------------------------------- +! Pi +pi=2.*asin(1.) + +!Reference surface pressure (Pa) +preff=101325. +CALL getin('preff', preff) +! Reference pressure at which hybrid coord. become purely pressure +! pa=50000. +pa=preff/2. +CALL getin('pa', pa) +! Gravity +g=9.80665 +CALL getin('g',g) +! Molar mass of the atmosphere +molmass = 28.9644 +CALL getin('molmass',molmass) +! kappa=R/Cp et Cp +kappa = 2./7. +CALL getin('kappa',kappa) +cpp=8.3145/molmass/kappa*1000. +CALL getin('cpp',cpp) +! Radius of the planet +rad = 6371229. +CALL getin('radius',rad) +! Length of a standard day (s) +daysec=86400. +CALL getin('daysec',daysec) +! Rotation rate of the planet: +! Length of a solar day, in standard days +daylen = 1. +CALL getin('daylen',daylen) +! Number of days (standard) per year: +year_day = 365.25 +CALL getin('year_day',year_day) +! Omega +! omeg=2.*pi/86400. +omeg=2.*pi/daysec*(1./daylen+1./year_day) +CALL getin('omeg',omeg) + +! Intrinsic heat flux (default: none) (only used if planet_type="giant") +ihf = 0. +call getin('ihf',ihf) + +END SUBROUTINE conf_planete Index: LMDZ5/trunk/libf/dyn3dmem/control_mod.F90 =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/control_mod.F90 (revision 1672) +++ LMDZ5/trunk/libf/dyn3dmem/control_mod.F90 (revision 1673) @@ -10,7 +10,7 @@ IMPLICIT NONE - REAL :: periodav + REAL :: periodav, starttime INTEGER :: nday,day_step,iperiod,iapp_tracvl,nsplit_phys - INTEGER :: iconser,iecri,idissip,iphysiq,iecrimoy + INTEGER :: iconser,iecri,dissip_period,iphysiq,iecrimoy INTEGER :: dayref,anneeref, raz_date, ip_ebil_dyn LOGICAL :: offline Index: LMDZ5/trunk/libf/dyn3dmem/defrun.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/defrun.F (revision 1672) +++ LMDZ5/trunk/libf/dyn3dmem/defrun.F (revision 1673) @@ -1,4 +1,4 @@ ! -! $Id: defrun.F 1403 2010-07-01 09:02:53Z fairhead $ +! $Id$ ! c @@ -132,7 +132,7 @@ READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) idissip - WRITE(tapeout,9001) ch1,'idissip' - WRITE(tapeout,*) idissip + READ (tapedef,*) dissip_period + WRITE(tapeout,9001) ch1,'dissip_period' + WRITE(tapeout,*) dissip_period ccc .... P. Le Van , modif le 29/04/97 .pour la dissipation ... Index: LMDZ5/trunk/libf/dyn3dmem/disvert.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/disvert.F (revision 1672) +++ (revision ) @@ -1,194 +1,0 @@ -! -! $Id: disvert.F 1403 2010-07-01 09:02:53Z fairhead $ -! - SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig) - -c Auteur : P. Le Van . -c - IMPLICIT NONE - -#include "dimensions.h" -#include "paramet.h" -#include "iniprint.h" -#include "logic.h" -c -c======================================================================= -c -c -c s = sigma ** kappa : coordonnee verticale -c dsig(l) : epaisseur de la couche l ds la coord. s -c sig(l) : sigma a l'interface des couches l et l-1 -c ds(l) : distance entre les couches l et l-1 en coord.s -c -c======================================================================= -c - REAL pa,preff - REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1) - REAL presnivs(llm) -c -c declarations: -c ------------- -c - REAL sig(llm+1),dsig(llm) - real zzz(1:llm+1) - real dzz(1:llm) - real zk,zkm1,dzk1,dzk2,k0,k1 -c - INTEGER l - REAL snorm,dsigmin - REAL alpha,beta,gama,delta,deltaz,h - INTEGER np,ierr - REAL pi,x - - REAL SSUM -c -c----------------------------------------------------------------------- -c - pi=2.*ASIN(1.) - - OPEN(99,file='sigma.def',status='old',form='formatted', - s iostat=ierr) - -c----------------------------------------------------------------------- -c cas 1 on lit les options dans sigma.def: -c ---------------------------------------- - - IF (ierr.eq.0) THEN - - READ(99,*) h ! hauteur d'echelle 8. - READ(99,*) deltaz ! epaiseur de la premiere couche 0.04 - READ(99,*) beta ! facteur d'acroissement en haut 1.3 - READ(99,*) k0 ! nombre de couches dans la transition surf - READ(99,*) k1 ! nombre de couches dans la transition haute - CLOSE(99) - alpha=deltaz/(llm*h) - write(lunout,*)'h,alpha,k0,k1,beta' - -c read(*,*) h,deltaz,beta,k0,k1 ! 8 0.04 4 20 1.2 - - alpha=deltaz/tanh(1./k0)*2. - zkm1=0. - sig(1)=1. - do l=1,llm - sig(l+1)=(cosh(l/k0))**(-alpha*k0/h) - + *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)) - zk=-h*log(sig(l+1)) - - dzk1=alpha*tanh(l/k0) - dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta) - write(lunout,*)l,sig(l+1),zk,zk-zkm1,dzk1,dzk2 - zkm1=zk - enddo - - sig(llm+1)=0. - -c - DO 2 l = 1, llm - dsig(l) = sig(l)-sig(l+1) - 2 CONTINUE -c - - ELSE -c----------------------------------------------------------------------- -c cas 2 ancienne discretisation (LMD5...): -c ---------------------------------------- - - WRITE(LUNOUT,*)'WARNING!!! Ancienne discretisation verticale' - - if (ok_strato) then - if (llm==39) then - dsigmin=0.3 - else if (llm==50) then - dsigmin=1. - else - WRITE(LUNOUT,*) 'ATTENTION discretisation z a ajuster' - dsigmin=1. - endif - WRITE(LUNOUT,*) 'Discretisation verticale DSIGMIN=',dsigmin - endif - - h=7. - snorm = 0. - DO l = 1, llm - x = 2.*asin(1.) * (REAL(l)-0.5) / REAL(llm+1) - - IF (ok_strato) THEN - dsig(l) =(dsigmin + 7.0 * SIN(x)**2) - & *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2 - ELSE - dsig(l) = 1.0 + 7.0 * SIN(x)**2 - ENDIF - - snorm = snorm + dsig(l) - ENDDO - snorm = 1./snorm - DO l = 1, llm - dsig(l) = dsig(l)*snorm - ENDDO - sig(llm+1) = 0. - DO l = llm, 1, -1 - sig(l) = sig(l+1) + dsig(l) - ENDDO - - ENDIF - - - DO l=1,llm - nivsigs(l) = REAL(l) - ENDDO - - DO l=1,llmp1 - nivsig(l)= REAL(l) - ENDDO - -c -c .... Calculs de ap(l) et de bp(l) .... -c ......................................... -c -c -c ..... pa et preff sont lus sur les fichiers start par lectba ..... -c - - bp(llmp1) = 0. - - DO l = 1, llm -cc -ccc ap(l) = 0. -ccc bp(l) = sig(l) - - bp(l) = EXP( 1. -1./( sig(l)*sig(l)) ) - ap(l) = pa * ( sig(l) - bp(l) ) -c - ENDDO - - bp(1)=1. - ap(1)=0. - - ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) ) - - write(lunout,*)' BP ' - write(lunout,*) bp - write(lunout,*)' AP ' - write(lunout,*) ap - - write(lunout,*) - .'Niveaux de pressions approximatifs aux centres des' - write(lunout,*)'couches calcules pour une pression de surface =', - . preff - write(lunout,*) - . 'et altitudes equivalentes pour une hauteur d echelle de' - write(lunout,*)'8km' - DO l = 1, llm - dpres(l) = bp(l) - bp(l+1) - presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff ) - write(lunout,*)'PRESNIVS(',l,')=',presnivs(l),' Z ~ ', - . log(preff/presnivs(l))*8. - . ,' DZ ~ ',8.*log((ap(l)+bp(l)*preff)/ - . max(ap(l+1)+bp(l+1)*preff,1.e-10)) - ENDDO - - write(lunout,*)' PRESNIVS ' - write(lunout,*)presnivs - - RETURN - END Index: LMDZ5/trunk/libf/dyn3dmem/disvert.F90 =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/disvert.F90 (revision 1673) +++ LMDZ5/trunk/libf/dyn3dmem/disvert.F90 (revision 1673) @@ -0,0 +1,180 @@ +! $Id: disvert.F90 1645 2012-07-30 16:01:50Z lguez $ + +SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,scaleheight) + + ! Auteur : P. Le Van + + use new_unit_m, only: new_unit + use ioipsl, only: getin + use assert_m, only: assert + + IMPLICIT NONE + + include "dimensions.h" + include "paramet.h" + include "iniprint.h" + include "logic.h" + + ! s = sigma ** kappa : coordonnee verticale + ! dsig(l) : epaisseur de la couche l ds la coord. s + ! sig(l) : sigma a l'interface des couches l et l-1 + ! ds(l) : distance entre les couches l et l-1 en coord.s + + real,intent(in) :: pa, preff + real,intent(out) :: ap(llmp1) ! in Pa + real, intent(out):: bp(llmp1) + real,intent(out) :: dpres(llm), nivsigs(llm), nivsig(llmp1) + real,intent(out) :: presnivs(llm) + real,intent(out) :: scaleheight + + REAL sig(llm+1), dsig(llm) + real zk, zkm1, dzk1, dzk2, k0, k1 + + INTEGER l, unit + REAL dsigmin + REAL alpha, beta, deltaz + REAL x + character(len=*),parameter :: modname="disvert" + + character(len=6):: vert_sampling + ! (allowed values are "param", "tropo", "strato" and "read") + + !----------------------------------------------------------------------- + + print *, "Call sequence information: disvert" + + ! default scaleheight is 8km for earth + scaleheight=8. + + vert_sampling = merge("strato", "tropo ", ok_strato) ! default value + call getin('vert_sampling', vert_sampling) + print *, 'vert_sampling = ' // vert_sampling + + select case (vert_sampling) + case ("param") + ! On lit les options dans sigma.def: + OPEN(99, file='sigma.def', status='old', form='formatted') + READ(99, *) scaleheight ! hauteur d'echelle 8. + READ(99, *) deltaz ! epaiseur de la premiere couche 0.04 + READ(99, *) beta ! facteur d'acroissement en haut 1.3 + READ(99, *) k0 ! nombre de couches dans la transition surf + READ(99, *) k1 ! nombre de couches dans la transition haute + CLOSE(99) + alpha=deltaz/(llm*scaleheight) + write(lunout, *)trim(modname),':scaleheight, alpha, k0, k1, beta', & + scaleheight, alpha, k0, k1, beta + + alpha=deltaz/tanh(1./k0)*2. + zkm1=0. + sig(1)=1. + do l=1, llm + sig(l+1)=(cosh(l/k0))**(-alpha*k0/scaleheight) & + *exp(-alpha/scaleheight*tanh((llm-k1)/k0) & + *beta**(l-(llm-k1))/log(beta)) + zk=-scaleheight*log(sig(l+1)) + + dzk1=alpha*tanh(l/k0) + dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta) + write(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2 + zkm1=zk + enddo + + sig(llm+1)=0. + + bp(: llm) = EXP(1. - 1. / sig(: llm)**2) + bp(llmp1) = 0. + + ap = pa * (sig - bp) + case("tropo") + DO l = 1, llm + x = 2*asin(1.) * (l - 0.5) / (llm + 1) + dsig(l) = 1.0 + 7.0 * SIN(x)**2 + ENDDO + dsig = dsig / sum(dsig) + sig(llm+1) = 0. + DO l = llm, 1, -1 + sig(l) = sig(l+1) + dsig(l) + ENDDO + + bp(1)=1. + bp(2: llm) = EXP(1. - 1. / sig(2: llm)**2) + bp(llmp1) = 0. + + ap(1)=0. + ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1)) + case("strato") + if (llm==39) then + dsigmin=0.3 + else if (llm==50) then + dsigmin=1. + else + write(lunout,*) trim(modname), ' ATTENTION discretisation z a ajuster' + dsigmin=1. + endif + WRITE(LUNOUT,*) trim(modname), 'Discretisation verticale DSIGMIN=',dsigmin + + DO l = 1, llm + x = 2*asin(1.) * (l - 0.5) / (llm + 1) + dsig(l) =(dsigmin + 7. * SIN(x)**2) & + *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2 + ENDDO + dsig = dsig / sum(dsig) + sig(llm+1) = 0. + DO l = llm, 1, -1 + sig(l) = sig(l+1) + dsig(l) + ENDDO + + bp(1)=1. + bp(2: llm) = EXP(1. - 1. / sig(2: llm)**2) + bp(llmp1) = 0. + + ap(1)=0. + ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1)) + case("read") + ! Read "ap" and "bp". First line is skipped (title line). "ap" + ! should be in Pa. First couple of values should correspond to + ! the surface, that is : "bp" should be in descending order. + call new_unit(unit) + open(unit, file="hybrid.txt", status="old", action="read", & + position="rewind") + read(unit, fmt=*) ! skip title line + do l = 1, llm + 1 + read(unit, fmt=*) ap(l), bp(l) + end do + close(unit) + call assert(ap(1) == 0., ap(llm + 1) == 0., bp(1) == 1., & + bp(llm + 1) == 0., "disvert: bad ap or bp values") + case default + call abort_gcm("disvert", 'Wrong value for "vert_sampling"', 1) + END select + + DO l=1, llm + nivsigs(l) = REAL(l) + ENDDO + + DO l=1, llmp1 + nivsig(l)= REAL(l) + ENDDO + + write(lunout, *) trim(modname),': BP ' + write(lunout, *) bp + write(lunout, *) trim(modname),': AP ' + write(lunout, *) ap + + write(lunout, *) 'Niveaux de pressions approximatifs aux centres des' + write(lunout, *)'couches calcules pour une pression de surface =', preff + write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de ' + write(lunout, *) scaleheight,' km' + DO l = 1, llm + dpres(l) = bp(l) - bp(l+1) + presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff ) + write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', & + log(preff/presnivs(l))*scaleheight & + , ' DZ ~ ', scaleheight*log((ap(l)+bp(l)*preff)/ & + max(ap(l+1)+bp(l+1)*preff, 1.e-10)) + ENDDO + + write(lunout, *) trim(modname),': PRESNIVS ' + write(lunout, *) presnivs + +END SUBROUTINE disvert Index: LMDZ5/trunk/libf/dyn3dmem/disvert_noterre.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/disvert_noterre.F (revision 1673) +++ LMDZ5/trunk/libf/dyn3dmem/disvert_noterre.F (revision 1673) @@ -0,0 +1,330 @@ +! $Id$ + SUBROUTINE disvert_noterre + +c Auteur : F. Forget Y. Wanherdrick, P. Levan +c Nouvelle version 100% Mars !! +c On l'utilise aussi pour Venus et Titan, legerment modifiee. + +#ifdef CPP_IOIPSL + use IOIPSL +#else +! if not using IOIPSL, we still need to use (a local version of) getin + use ioipsl_getincom +#endif + + IMPLICIT NONE + +#include "dimensions.h" +#include "paramet.h" +#include "comvert.h" +#include "comconst.h" +#include "logic.h" +#include "iniprint.h" +c +c======================================================================= +c Discretisation verticale en coordonnée hybride (ou sigma) +c +c======================================================================= +c +c declarations: +c ------------- +c +c + INTEGER l,ll + REAL snorm + REAL alpha,beta,gama,delta,deltaz + real quoi,quand + REAL zsig(llm),sig(llm+1) + INTEGER np,ierr + integer :: ierr1,ierr2,ierr3,ierr4 + REAL x + + REAL SSUM + EXTERNAL SSUM + real newsig + REAL dz0,dz1,nhaut,sig1,esig,csig,zz + real tt,rr,gg, prevz + real s(llm),dsig(llm) + + integer iz + real z, ps,p + character(len=*),parameter :: modname="disvert_noterre" + +c +c----------------------------------------------------------------------- +c +! Initializations: +! pi=2.*ASIN(1.) ! already done in iniconst + + hybrid=.true. ! default value for hybrid (ie: use hybrid coordinates) + CALL getin('hybrid',hybrid) + write(lunout,*) trim(modname),': hybrid=',hybrid + +! Ouverture possible de fichiers typiquement E.T. + + open(99,file="esasig.def",status='old',form='formatted', + s iostat=ierr2) + if(ierr2.ne.0) then + close(99) + open(99,file="z2sig.def",status='old',form='formatted', + s iostat=ierr4) + endif + +c----------------------------------------------------------------------- +c cas 1 on lit les options dans esasig.def: +c ---------------------------------------- + + IF(ierr2.eq.0) then + +c Lecture de esasig.def : +c Systeme peu souple, mais qui respecte en theorie +c La conservation de l'energie (conversion Energie potentielle +c <-> energie cinetique, d'apres la note de Frederic Hourdin... + + write(lunout,*)'*****************************' + write(lunout,*)'WARNING reading esasig.def' + write(lunout,*)'*****************************' + READ(99,*) scaleheight + READ(99,*) dz0 + READ(99,*) dz1 + READ(99,*) nhaut + CLOSE(99) + + dz0=dz0/scaleheight + dz1=dz1/scaleheight + + sig1=(1.-dz1)/tanh(.5*(llm-1)/nhaut) + + esig=1. + + do l=1,20 + esig=-log((1./sig1-1.)*exp(-dz0)/esig)/(llm-1.) + enddo + csig=(1./sig1-1.)/(exp(esig)-1.) + + DO L = 2, llm + zz=csig*(exp(esig*(l-1.))-1.) + sig(l) =1./(1.+zz) + & * tanh(.5*(llm+1-l)/nhaut) + ENDDO + sig(1)=1. + sig(llm+1)=0. + quoi = 1. + 2.* kappa + s( llm ) = 1. + s(llm-1) = quoi + IF( llm.gt.2 ) THEN + DO ll = 2, llm-1 + l = llm+1 - ll + quand = sig(l+1)/ sig(l) + s(l-1) = quoi * (1.-quand) * s(l) + quand * s(l+1) + ENDDO + END IF +c + snorm=(1.-.5*sig(2)+kappa*(1.-sig(2)))*s(1)+.5*sig(2)*s(2) + DO l = 1, llm + s(l) = s(l)/ snorm + ENDDO + +c----------------------------------------------------------------------- +c cas 2 on lit les options dans z2sig.def: +c ---------------------------------------- + + ELSE IF(ierr4.eq.0) then + write(lunout,*)'****************************' + write(lunout,*)'Reading z2sig.def' + write(lunout,*)'****************************' + + READ(99,*) scaleheight + do l=1,llm + read(99,*) zsig(l) + end do + CLOSE(99) + + sig(1) =1 + do l=2,llm + sig(l) = 0.5 * ( exp(-zsig(l)/scaleheight) + + & exp(-zsig(l-1)/scaleheight) ) + end do + sig(llm+1) =0 + +c----------------------------------------------------------------------- + ELSE + write(lunout,*) 'didn t you forget something ??? ' + write(lunout,*) 'We need file z2sig.def ! (OR esasig.def)' + stop + ENDIF +c----------------------------------------------------------------------- + + DO l=1,llm + nivsigs(l) = REAL(l) + ENDDO + + DO l=1,llmp1 + nivsig(l)= REAL(l) + ENDDO + + +c----------------------------------------------------------------------- +c .... Calculs de ap(l) et de bp(l) .... +c ......................................... +c +c ..... pa et preff sont lus sur les fichiers start par dynetat0 ..... +c----------------------------------------------------------------------- +c + + if (hybrid) then ! use hybrid coordinates + write(lunout,*) "*********************************" + write(lunout,*) "Using hybrid vertical coordinates" + write(lunout,*) +c Coordonnees hybrides avec mod + DO l = 1, llm + + call sig_hybrid(sig(l),pa,preff,newsig) + bp(l) = EXP( 1. - 1./(newsig**2) ) + ap(l) = pa * (newsig - bp(l) ) + enddo + bp(llmp1) = 0. + ap(llmp1) = 0. + else ! use sigma coordinates + write(lunout,*) "********************************" + write(lunout,*) "Using sigma vertical coordinates" + write(lunout,*) +c Pour ne pas passer en coordonnees hybrides + DO l = 1, llm + ap(l) = 0. + bp(l) = sig(l) + ENDDO + ap(llmp1) = 0. + endif + + bp(llmp1) = 0. + + write(lunout,*) trim(modname),': BP ' + write(lunout,*) bp + write(lunout,*) trim(modname),': AP ' + write(lunout,*) ap + +c Calcul au milieu des couches : +c WARNING : le choix de placer le milieu des couches au niveau de +c pression intermédiaire est arbitraire et pourrait etre modifié. +c Le calcul du niveau pour la derniere couche +c (on met la meme distance (en log pression) entre P(llm) +c et P(llm -1) qu'entre P(llm-1) et P(llm-2) ) est +c Specifique. Ce choix est spécifié ici ET dans exner_milieu.F + + DO l = 1, llm-1 + aps(l) = 0.5 *( ap(l) +ap(l+1)) + bps(l) = 0.5 *( bp(l) +bp(l+1)) + ENDDO + + if (hybrid) then + aps(llm) = aps(llm-1)**2 / aps(llm-2) + bps(llm) = 0.5*(bp(llm) + bp(llm+1)) + else + bps(llm) = bps(llm-1)**2 / bps(llm-2) + aps(llm) = 0. + end if + + write(lunout,*) trim(modname),': BPs ' + write(lunout,*) bps + write(lunout,*) trim(modname),': APs' + write(lunout,*) aps + + DO l = 1, llm + presnivs(l) = aps(l)+bps(l)*preff + pseudoalt(l) = -scaleheight*log(presnivs(l)/preff) + ENDDO + + write(lunout,*)trim(modname),' : PRESNIVS' + write(lunout,*)presnivs + write(lunout,*)'Pseudo altitude of Presnivs : (for a scale ', + & 'height of ',scaleheight,' km)' + write(lunout,*)pseudoalt + +c -------------------------------------------------- +c This can be used to plot the vertical discretization +c (> xmgrace -nxy testhybrid.tab ) +c -------------------------------------------------- +c open (53,file='testhybrid.tab') +c scaleheight=15.5 +c do iz=0,34 +c z = -5 + min(iz,34-iz) +c approximation of scale height for Venus +c scaleheight = 15.5 - z/55.*10. +c ps = preff*exp(-z/scaleheight) +c zsig(1)= -scaleheight*log((aps(1) + bps(1)*ps)/preff) +c do l=2,llm +c approximation of scale height for Venus +c if (zsig(l-1).le.55.) then +c scaleheight = 15.5 - zsig(l-1)/55.*10. +c else +c scaleheight = 5.5 - (zsig(l-1)-55.)/35.*2. +c endif +c zsig(l)= zsig(l-1)-scaleheight* +c . log((aps(l) + bps(l)*ps)/(aps(l-1) + bps(l-1)*ps)) +c end do +c write(53,'(I3,50F10.5)') iz, zsig +c end do +c close(53) +c -------------------------------------------------- + + + RETURN + END + +c ************************************************************ + subroutine sig_hybrid(sig,pa,preff,newsig) +c ---------------------------------------------- +c Subroutine utilisee pour calculer des valeurs de sigma modifie +c pour conserver les coordonnees verticales decrites dans +c esasig.def/z2sig.def lors du passage en coordonnees hybrides +c F. Forget 2002 +c Connaissant sig (niveaux "sigma" ou on veut mettre les couches) +c L'objectif est de calculer newsig telle que +c (1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig = sig +c Cela ne se résoud pas analytiquement: +c => on résoud par iterration bourrine +c ---------------------------------------------- +c Information : where exp(1-1./x**2) become << x +c x exp(1-1./x**2) /x +c 1 1 +c 0.68 0.5 +c 0.5 1.E-1 +c 0.391 1.E-2 +c 0.333 1.E-3 +c 0.295 1.E-4 +c 0.269 1.E-5 +c 0.248 1.E-6 +c => on peut utiliser newsig = sig*preff/pa si sig*preff/pa < 0.25 + + + implicit none + real x1, x2, sig,pa,preff, newsig, F + integer j + + newsig = sig + x1=0 + x2=1 + if (sig.ge.1) then + newsig= sig + else if (sig*preff/pa.ge.0.25) then + DO J=1,9999 ! nombre d''iteration max + F=((1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig)/sig +c write(0,*) J, ' newsig =', newsig, ' F= ', F + if (F.gt.1) then + X2 = newsig + newsig=(X1+newsig)*0.5 + else + X1 = newsig + newsig=(X2+newsig)*0.5 + end if +c Test : on arete lorsque on approxime sig à moins de 0.01 m près +c (en pseudo altitude) : + IF(abs(10.*log(F)).LT.1.E-5) goto 999 + END DO + else ! if (sig*preff/pa.le.0.25) then + newsig= sig*preff/pa + end if + 999 continue + Return + END Index: LMDZ5/trunk/libf/dyn3dmem/dynetat0.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/dynetat0.F (revision 1672) +++ (revision ) @@ -1,386 +1,0 @@ -! -! $Id $ -! - SUBROUTINE dynetat0(fichnom,vcov,ucov, - . teta,q,masse,ps,phis,time) - - USE infotrac - IMPLICIT NONE - -c======================================================================= -c -c Auteur: P. Le Van / L.Fairhead -c ------- -c -c objet: -c ------ -c -c Lecture de l'etat initial -c -c======================================================================= -c----------------------------------------------------------------------- -c Declarations: -c ------------- - -#include "dimensions.h" -#include "paramet.h" -#include "temps.h" -#include "comconst.h" -#include "comvert.h" -#include "comgeom.h" -#include "ener.h" -#include "netcdf.inc" -#include "description.h" -#include "serre.h" -#include "logic.h" -#include "iniprint.h" - -c Arguments: -c ---------- - - CHARACTER*(*) fichnom - REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm),teta(ip1jmp1,llm) - REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm) - REAL ps(ip1jmp1),phis(ip1jmp1) - - REAL time - -c Variables -c - INTEGER length,iq - PARAMETER (length = 100) - REAL tab_cntrl(length) ! tableau des parametres du run - INTEGER ierr, nid, nvarid - -c----------------------------------------------------------------------- - -c Ouverture NetCDF du fichier etat initial - - ierr = NF_OPEN (fichnom, NF_NOWRITE,nid) - IF (ierr.NE.NF_NOERR) THEN - write(lunout,*)'dynetat0: Pb d''ouverture du fichier start.nc' - write(lunout,*)' ierr = ', ierr - CALL ABORT - ENDIF - -c - ierr = NF_INQ_VARID (nid, "controle", nvarid) - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Le champ est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl) -#else - ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl) -#endif - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Lecture echoue pour " - CALL abort - ENDIF - - im = tab_cntrl(1) - jm = tab_cntrl(2) - lllm = tab_cntrl(3) - day_ref = tab_cntrl(4) - annee_ref = tab_cntrl(5) - rad = tab_cntrl(6) - omeg = tab_cntrl(7) - g = tab_cntrl(8) - cpp = tab_cntrl(9) - kappa = tab_cntrl(10) - daysec = tab_cntrl(11) - dtvr = tab_cntrl(12) - etot0 = tab_cntrl(13) - ptot0 = tab_cntrl(14) - ztot0 = tab_cntrl(15) - stot0 = tab_cntrl(16) - ang0 = tab_cntrl(17) - pa = tab_cntrl(18) - preff = tab_cntrl(19) -c - clon = tab_cntrl(20) - clat = tab_cntrl(21) - grossismx = tab_cntrl(22) - grossismy = tab_cntrl(23) -c - IF ( tab_cntrl(24).EQ.1. ) THEN - fxyhypb = . TRUE . -c dzoomx = tab_cntrl(25) -c dzoomy = tab_cntrl(26) -c taux = tab_cntrl(28) -c tauy = tab_cntrl(29) - ELSE - fxyhypb = . FALSE . - ysinus = . FALSE . - IF( tab_cntrl(27).EQ.1. ) ysinus = . TRUE. - ENDIF - - day_ini = tab_cntrl(30) - itau_dyn = tab_cntrl(31) -c ................................................................. -c -c - write(lunout,*)'dynetat0: rad,omeg,g,cpp,kappa', - & rad,omeg,g,cpp,kappa - - IF( im.ne.iim ) THEN - PRINT 1,im,iim - STOP - ELSE IF( jm.ne.jjm ) THEN - PRINT 2,jm,jjm - STOP - ELSE IF( lllm.ne.llm ) THEN - PRINT 3,lllm,llm - STOP - ENDIF - - ierr = NF_INQ_VARID (nid, "rlonu", nvarid) - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Le champ est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonu) -#else - ierr = NF_GET_VAR_REAL(nid, nvarid, rlonu) -#endif - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Lecture echouee pour " - CALL abort - ENDIF - - ierr = NF_INQ_VARID (nid, "rlatu", nvarid) - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Le champ est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatu) -#else - ierr = NF_GET_VAR_REAL(nid, nvarid, rlatu) -#endif - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Lecture echouee pour " - CALL abort - ENDIF - - ierr = NF_INQ_VARID (nid, "rlonv", nvarid) - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Le champ est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonv) -#else - ierr = NF_GET_VAR_REAL(nid, nvarid, rlonv) -#endif - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Lecture echouee pour " - CALL abort - ENDIF - - ierr = NF_INQ_VARID (nid, "rlatv", nvarid) - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Le champ est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatv) -#else - ierr = NF_GET_VAR_REAL(nid, nvarid, rlatv) -#endif - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Lecture echouee pour rlatv" - CALL abort - ENDIF - - ierr = NF_INQ_VARID (nid, "cu", nvarid) - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Le champ est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VAR_DOUBLE(nid, nvarid, cu) -#else - ierr = NF_GET_VAR_REAL(nid, nvarid, cu) -#endif - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Lecture echouee pour " - CALL abort - ENDIF - - ierr = NF_INQ_VARID (nid, "cv", nvarid) - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Le champ est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VAR_DOUBLE(nid, nvarid, cv) -#else - ierr = NF_GET_VAR_REAL(nid, nvarid, cv) -#endif - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Lecture echouee pour " - CALL abort - ENDIF - - ierr = NF_INQ_VARID (nid, "aire", nvarid) - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Le champ est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VAR_DOUBLE(nid, nvarid, aire) -#else - ierr = NF_GET_VAR_REAL(nid, nvarid, aire) -#endif - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Lecture echouee pour " - CALL abort - ENDIF - - ierr = NF_INQ_VARID (nid, "phisinit", nvarid) - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Le champ est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VAR_DOUBLE(nid, nvarid, phis) -#else - ierr = NF_GET_VAR_REAL(nid, nvarid, phis) -#endif - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Lecture echouee pour " - CALL abort - ENDIF - - ierr = NF_INQ_VARID (nid, "temps", nvarid) - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Le champ est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VAR_DOUBLE(nid, nvarid, time) -#else - ierr = NF_GET_VAR_REAL(nid, nvarid, time) -#endif - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Lecture echouee " - CALL abort - ENDIF - - ierr = NF_INQ_VARID (nid, "ucov", nvarid) - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Le champ est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ucov) -#else - ierr = NF_GET_VAR_REAL(nid, nvarid, ucov) -#endif - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Lecture echouee pour " - CALL abort - ENDIF - - ierr = NF_INQ_VARID (nid, "vcov", nvarid) - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Le champ est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VAR_DOUBLE(nid, nvarid, vcov) -#else - ierr = NF_GET_VAR_REAL(nid, nvarid, vcov) -#endif - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Lecture echouee pour " - CALL abort - ENDIF - - ierr = NF_INQ_VARID (nid, "teta", nvarid) - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Le champ est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VAR_DOUBLE(nid, nvarid, teta) -#else - ierr = NF_GET_VAR_REAL(nid, nvarid, teta) -#endif - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Lecture echouee pour " - CALL abort - ENDIF - - - IF(nqtot.GE.1) THEN - DO iq=1,nqtot - ierr = NF_INQ_VARID (nid, tname(iq), nvarid) - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Le traceur <"//trim(tname(iq))// - & "> est absent" - write(lunout,*)" Il est donc initialise a zero" - q(:,:,iq)=0. - ELSE -#ifdef NC_DOUBLE - ierr = NF_GET_VAR_DOUBLE(nid, nvarid, q(1,1,iq)) -#else - ierr = NF_GET_VAR_REAL(nid, nvarid, q(1,1,iq)) -#endif - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Lecture echouee pour "//tname(iq) - CALL abort - ENDIF - ENDIF - ENDDO - ENDIF - - ierr = NF_INQ_VARID (nid, "masse", nvarid) - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Le champ est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VAR_DOUBLE(nid, nvarid, masse) -#else - ierr = NF_GET_VAR_REAL(nid, nvarid, masse) -#endif - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Lecture echouee pour " - CALL abort - ENDIF - - ierr = NF_INQ_VARID (nid, "ps", nvarid) - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Le champ est absent" - CALL abort - ENDIF -#ifdef NC_DOUBLE - ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ps) -#else - ierr = NF_GET_VAR_REAL(nid, nvarid, ps) -#endif - IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0: Lecture echouee pour " - CALL abort - ENDIF - - ierr = NF_CLOSE(nid) - - day_ini=day_ini+INT(time) - time=time-INT(time) - - 1 FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dem - *arrage est differente de la valeur parametree iim =',i4//) - 2 FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dem - *arrage est differente de la valeur parametree jjm =',i4//) - 3 FORMAT(//10x,'la valeur de lmax =',i4,2x,'lue sur le fichier dema - *rrage est differente de la valeur parametree llm =',i4//) - 4 FORMAT(//10x,'la valeur de dtrv =',i4,2x,'lue sur le fichier dema - *rrage est differente de la valeur dtinteg =',i4//) - - RETURN - END Index: LMDZ5/trunk/libf/dyn3dmem/dynetat0_loc.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/dynetat0_loc.F (revision 1672) +++ LMDZ5/trunk/libf/dyn3dmem/dynetat0_loc.F (revision 1673) @@ -5,4 +5,5 @@ . teta,q,masse,ps,phis,time) USE infotrac + use control_mod, only : planet_type USE parallel IMPLICIT NONE @@ -57,4 +58,6 @@ REAL,ALLOCATABLE :: phis_glo(:) + INTEGER idecal + c----------------------------------------------------------------------- c Ouverture NetCDF du fichier etat initial @@ -84,30 +87,42 @@ ENDIF + !!! AS: idecal is a hack to be able to read planeto starts... + !!! .... while keeping everything OK for LMDZ EARTH + if (planet_type.eq."generic") then + print*,'NOTE NOTE NOTE : Planeto-like start files' + idecal = 4 + annee_ref = 2000 + else + print*,'NOTE NOTE NOTE : Earth-like start files' + idecal = 5 + annee_ref = tab_cntrl(5) + endif + + im = tab_cntrl(1) jm = tab_cntrl(2) lllm = tab_cntrl(3) day_ref = tab_cntrl(4) - annee_ref = tab_cntrl(5) - rad = tab_cntrl(6) - omeg = tab_cntrl(7) - g = tab_cntrl(8) - cpp = tab_cntrl(9) - kappa = tab_cntrl(10) - daysec = tab_cntrl(11) - dtvr = tab_cntrl(12) - etot0 = tab_cntrl(13) - ptot0 = tab_cntrl(14) - ztot0 = tab_cntrl(15) - stot0 = tab_cntrl(16) - ang0 = tab_cntrl(17) - pa = tab_cntrl(18) - preff = tab_cntrl(19) -c - clon = tab_cntrl(20) - clat = tab_cntrl(21) - grossismx = tab_cntrl(22) - grossismy = tab_cntrl(23) -c - IF ( tab_cntrl(24).EQ.1. ) THEN + rad = tab_cntrl(idecal+1) + omeg = tab_cntrl(idecal+2) + g = tab_cntrl(idecal+3) + cpp = tab_cntrl(idecal+4) + kappa = tab_cntrl(idecal+5) + daysec = tab_cntrl(idecal+6) + dtvr = tab_cntrl(idecal+7) + etot0 = tab_cntrl(idecal+8) + ptot0 = tab_cntrl(idecal+9) + ztot0 = tab_cntrl(idecal+10) + stot0 = tab_cntrl(idecal+11) + ang0 = tab_cntrl(idecal+12) + pa = tab_cntrl(idecal+13) + preff = tab_cntrl(idecal+14) +c + clon = tab_cntrl(idecal+15) + clat = tab_cntrl(idecal+16) + grossismx = tab_cntrl(idecal+17) + grossismy = tab_cntrl(idecal+18) +c + IF ( tab_cntrl(idecal+19).EQ.1. ) THEN fxyhypb = . TRUE . c dzoomx = tab_cntrl(25) @@ -118,5 +133,5 @@ fxyhypb = . FALSE . ysinus = . FALSE . - IF( tab_cntrl(27).EQ.1. ) ysinus = . TRUE. + IF( tab_cntrl(idecal+22).EQ.1. ) ysinus = . TRUE. ENDIF @@ -266,6 +281,11 @@ ierr = NF_INQ_VARID (nid, "temps", nvarid) IF (ierr .NE. NF_NOERR) THEN - write(lunout,*)"dynetat0_loc: Le champ est absent" - CALL abort + write(lunout,*)"dynetat0: Le champ est absent" + write(lunout,*)"dynetat0: J essaie