Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3d/comconst.h =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/dyn3d/comconst.h (revision 1189) +++ LMDZ4/branches/LMDZ4-dev/libf/dyn3d/comconst.h (revision 1190) @@ -8,5 +8,6 @@ & dtvr,daysec, & & pi,dtphys,dtdiss,rad,r,cpp,kappa,cotot,unsim,g,omeg & - & ,dissip_factz,dissip_deltaz,dissip_zref + & ,dissip_factz,dissip_deltaz,dissip_zref & + & ,iflag_top_bound,tau_top_bound @@ -16,4 +17,7 @@ REAL cotot,unsim,g,omeg REAL dissip_factz,dissip_deltaz,dissip_zref + INTEGER iflag_top_bound + REAL tau_top_bound + !----------------------------------------------------------------------- Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3d/conf_gcm.F =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/dyn3d/conf_gcm.F (revision 1189) +++ LMDZ4/branches/LMDZ4-dev/libf/dyn3d/conf_gcm.F (revision 1190) @@ -280,4 +280,9 @@ CALL getin('dissip_deltaz',dissip_deltaz ) CALL getin('dissip_zref',dissip_zref ) + + iflag_top_bound=1 + tau_top_bound=1.e-5 + CALL getin('iflag_top_bound',iflag_top_bound) + CALL getin('tau_top_bound',tau_top_bound) !Config Key = coefdis Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3d/leapfrog.F =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/dyn3d/leapfrog.F (revision 1189) +++ LMDZ4/branches/LMDZ4-dev/libf/dyn3d/leapfrog.F (revision 1190) @@ -386,5 +386,5 @@ IF (ok_strato) THEN - CALL top_bound( vcov,ucov,teta, dufi,dvfi,dtetafi) + CALL top_bound( vcov,ucov,teta,masse,dufi,dvfi,dtetafi) ENDIF Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3d/top_bound.F =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/dyn3d/top_bound.F (revision 1189) +++ LMDZ4/branches/LMDZ4-dev/libf/dyn3d/top_bound.F (revision 1190) @@ -1,3 +1,3 @@ - SUBROUTINE top_bound( vcov,ucov,teta, du,dv,dh ) + SUBROUTINE top_bound( vcov,ucov,teta,masse, du,dv,dh ) IMPLICIT NONE c @@ -5,5 +5,6 @@ #include "paramet.h" #include "comconst.h" -CC#include "comgeom2.h" +#include "comvert.h" +#include "comgeom2.h" @@ -27,5 +28,5 @@ c ------------- -#include "comgeom.h" +! #include "comgeom.h" #include "comdissipn.h" @@ -34,4 +35,5 @@ REAL ucov(iip1,jjp1,llm),vcov(iip1,jjm,llm),teta(iip1,jjp1,llm) + REAL masse(iip1,jjp1,llm) REAL dv(iip1,jjm,llm),du(iip1,jjp1,llm),dh(iip1,jjp1,llm) @@ -39,11 +41,14 @@ c ------ + REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm),zm REAL uzon(jjp1,llm),vzon(jjm,llm),tzon(jjp1,llm) INTEGER NDAMP PARAMETER (NDAMP=4) - integer i - REAL :: rdamp(llm) = - & (/(0., i =1,llm-NDAMP),0.125E-5,.25E-5,.5E-5,1.E-5/) + integer i + REAL,SAVE :: rdamp(llm) +! & (/(0., i =1,llm-NDAMP),0.125E-5,.25E-5,.5E-5,1.E-5/) + + LOGICAL,SAVE :: first=.true. INTEGER j,l @@ -52,10 +57,38 @@ C CALCUL DES CHAMPS EN MOYENNE ZONALE: + if (iflag_top_bound.eq.0) return + + if (first) then + if (iflag_top_bound.eq.1) then +! couche eponge dans les 4 dernieres couches du modele + rdamp(:)=0. + rdamp(llm)=tau_top_bound + rdamp(llm-1)=tau_top_bound/2. + rdamp(llm-2)=tau_top_bound/4. + rdamp(llm-3)=tau_top_bound/8. + else if (iflag_top_bound.eq.2) then +! couce eponge dans toutes les couches de pression plus faible que +! 100 fois la pression de la derniere couche + rdamp(:)=tau_top_bound + s *max(presnivs(llm)/presnivs(:)-0.01,0.) + endif + first=.false. + print*,'TOP_BOUND rdamp=',rdamp + endif + + CALL massbar(masse,massebx,masseby) + do l=1,llm do j=1,jjm vzon(j,l)=0. + zm=0. do i=1,iim - vzon(j,l)=vzon(j,l)+vcov(i,j,l)/float(iim) +! Rm: on peut travailler directement avec la moyenne zonale de vcov +! plutot qu'avec celle de v car le coefficient cv qui relie les deux +! ne varie qu'en latitude + vzon(j,l)=vzon(j,l)+vcov(i,j,l)*masseby(i,j,l) + zm=zm+masseby(i,j,l) enddo + vzon(j,l)=vzon(j,l)/zm enddo enddo @@ -72,9 +105,22 @@ do j=2,jjm uzon(j,l)=0. + zm=0. + do i=1,iim + uzon(j,l)=uzon(j,l)+massebx(i,j,l)*ucov(i,j,l)/cu(i,j) + zm=zm+massebx(i,j,l) + enddo + uzon(j,l)=uzon(j,l)/zm + enddo + enddo + + do l=1,llm + do j=2,jjm + zm=0. tzon(j,l)=0. do i=1,iim - uzon(j,l)=uzon(j,l)+ucov(i,j,l)/float(iim) - tzon(j,l)=tzon(j,l)+teta(i,j,l)/float(iim) + tzon(j,l)=tzon(j,l)+teta(i,j,l)*masse(i,j,l) + zm=zm+masse(i,j,l) enddo + tzon(j,l)=tzon(j,l)/zm enddo enddo @@ -85,5 +131,6 @@ do i=1,iip1 do j=2,jjm - du(i,j,l)=du(i,j,l)-rdamp(l)*(ucov(i,j,l)-uzon(j,l)) + du(i,j,l)=du(i,j,l) + s -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l)) dh(i,j,l)=dh(i,j,l)-rdamp(l)*(teta(i,j,l)-tzon(j,l)) enddo Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/comconst.h =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/comconst.h (revision 1189) +++ LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/comconst.h (revision 1190) @@ -8,5 +8,6 @@ & dtvr,daysec, & & pi,dtphys,dtdiss,rad,r,cpp,kappa,cotot,unsim,g,omeg & - & ,dissip_factz,dissip_deltaz,dissip_zref + & ,dissip_factz,dissip_deltaz,dissip_zref & + & ,iflag_top_bound,tau_top_bound @@ -16,4 +17,7 @@ REAL cotot,unsim,g,omeg REAL dissip_factz,dissip_deltaz,dissip_zref + INTEGER iflag_top_bound + REAL tau_top_bound + !----------------------------------------------------------------------- Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/conf_gcm.F =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/conf_gcm.F (revision 1189) +++ LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/conf_gcm.F (revision 1190) @@ -291,4 +291,10 @@ CALL getin('dissip_zref',dissip_zref ) + iflag_top_bound=1 + tau_top_bound=1.e-5 + CALL getin('iflag_top_bound',iflag_top_bound) + CALL getin('tau_top_bound',tau_top_bound) + +! !Config Key = coefdis !Config Desc = coefficient pour gamdissip Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/leapfrog_p.F =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/leapfrog_p.F (revision 1189) +++ LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/leapfrog_p.F (revision 1190) @@ -722,4 +722,7 @@ * jj_Nb_physic,2,2,Request_physic) + call Register_SwapFieldHallo(masse,masse,ip1jmp1,llm, + * jj_Nb_physic,1,2,Request_physic) + call Register_SwapFieldHallo(p,p,ip1jmp1,llmp1, * jj_Nb_physic,2,2,Request_physic) @@ -863,5 +866,5 @@ c ------------------------------ IF (ok_strato) THEN - CALL top_bound_p( vcov,ucov,teta, dufi,dvfi,dtetafi) + CALL top_bound_p( vcov,ucov,teta,masse,dufi,dvfi,dtetafi) ENDIF @@ -887,4 +890,7 @@ * jj_Nb_caldyn,Request_physic) + call Register_SwapField(masse,masse,ip1jmp1,llm, + * jj_Nb_caldyn,Request_physic) + call Register_SwapField(p,p,ip1jmp1,llmp1, * jj_Nb_caldyn,Request_physic) Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/top_bound_p.F =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/top_bound_p.F (revision 1189) +++ LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/top_bound_p.F (revision 1190) @@ -1,3 +1,3 @@ - SUBROUTINE top_bound_p( vcov,ucov,teta, du,dv,dh ) + SUBROUTINE top_bound_p( vcov,ucov,teta,masse, du,dv,dh ) USE parallel IMPLICIT NONE @@ -6,5 +6,6 @@ #include "paramet.h" #include "comconst.h" -CC#include "comgeom2.h" +#include "comvert.h" +#include "comgeom2.h" @@ -28,5 +29,4 @@ c ------------- -#include "comgeom.h" #include "comdissipn.h" @@ -35,9 +35,10 @@ REAL ucov(iip1,jjp1,llm),vcov(iip1,jjm,llm),teta(iip1,jjp1,llm) + REAL masse(iip1,jjp1,llm) REAL dv(iip1,jjm,llm),du(iip1,jjp1,llm),dh(iip1,jjp1,llm) c Local: c ------ - + REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm),zm REAL uzon(jjp1,llm),vzon(jjm,llm),tzon(jjp1,llm) @@ -45,10 +46,35 @@ PARAMETER (NDAMP=4) integer i - REAL :: rdamp(llm) = - & (/(0., i =1,llm-NDAMP),0.125E-5,.25E-5,.5E-5,1.E-5/) - + REAL,SAVE :: rdamp(llm) +! & (/(0., i =1,llm-NDAMP),0.125E-5,.25E-5,.5E-5,1.E-5/) + LOGICAL,SAVE :: first=.true. INTEGER j,l,jjb,jje + if (iflag_top_bound == 0) return + if (first) then +c$OMP BARRIER +c$OMP MASTER + if (iflag_top_bound == 1) then +! couche eponge dans les 4 dernieres couches du modele + rdamp(:)=0. + rdamp(llm)=tau_top_bound + rdamp(llm-1)=tau_top_bound/2. + rdamp(llm-2)=tau_top_bound/4. + rdamp(llm-3)=tau_top_bound/8. + else if (iflag_top_bound == 2) then +! couce eponge dans toutes les couches de pression plus faible que +! 100 fois la pression de la derniere couche + rdamp(:)=tau_top_bound + s *max(presnivs(llm)/presnivs(:)-0.01,0.) + endif + first=.false. + print*,'TOP_BOUND rdamp=',rdamp +c$OMP END MASTER +c$OMP BARRIER + endif + + + CALL massbar_p(masse,massebx,masseby) C CALCUL DES CHAMPS EN MOYENNE ZONALE: @@ -60,8 +86,14 @@ do l=1,llm do j=jjb,jje - vzon(j,l)=0. + zm=0. + vzon(j,l)=0 do i=1,iim - vzon(j,l)=vzon(j,l)+vcov(i,j,l)/float(iim) +! Rm: on peut travailler directement avec la moyenne zonale de vcov +! plutot qu'avec celle de v car le coefficient cv qui relie les deux +! ne varie qu'en latitude + vzon(j,l)=vzon(j,l)+vcov(i,j,l)*masseby(i,j,l) + zm=zm+masseby(i,j,l) enddo + vzon(j,l)=vzon(j,l)/zm enddo enddo @@ -87,9 +119,24 @@ do j=jjb,jje uzon(j,l)=0. + zm=0. + do i=1,iim + uzon(j,l)=uzon(j,l)+massebx(i,j,l)*ucov(i,j,l)/cu(i,j) + zm=zm+massebx(i,j,l) + enddo + uzon(j,l)=uzon(j,l)/zm + enddo + enddo +c$OMP END DO NOWAIT + +c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) + do l=1,llm + do j=jjb,jje + zm=0. tzon(j,l)=0. do i=1,iim - uzon(j,l)=uzon(j,l)+ucov(i,j,l)/float(iim) - tzon(j,l)=tzon(j,l)+teta(i,j,l)/float(iim) + tzon(j,l)=tzon(j,l)+teta(i,j,l)*masse(i,j,l) + zm=zm+masse(i,j,l) enddo + tzon(j,l)=tzon(j,l)/zm enddo enddo @@ -102,5 +149,6 @@ do j=jjb,jje do i=1,iip1 - du(i,j,l)=du(i,j,l)-rdamp(l)*(ucov(i,j,l)-uzon(j,l)) + du(i,j,l)=du(i,j,l) + s -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l)) dh(i,j,l)=dh(i,j,l)-rdamp(l)*(teta(i,j,l)-tzon(j,l)) enddo