Changeset 1190 for LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar

Timestamp:

Jun 23, 2009, 11:35:20 PM (15 years ago)

Author:

yann meurdesoif

Message:

Correction de FH sur la sponge layer de top_bound + parallelisation des corrections.

YM

Location:

LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar

Files:

• comconst.h (modified) (2 diffs)
• conf_gcm.F (modified) (1 diff)
• leapfrog_p.F (modified) (3 diffs)
• top_bound_p.F (modified) (8 diffs)

LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/comconst.h

@@ -8 +8 @@
      & 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 +17 @@
       REAL cotot,unsim,g,omeg
       REAL dissip_factz,dissip_deltaz,dissip_zref
+      INTEGER iflag_top_bound
+      REAL tau_top_bound
+
 
 !-----------------------------------------------------------------------

LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/conf_gcm.F

@@ -291 +291 @@
        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

LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/leapfrog_p.F

@@ -722 +722 @@
      *                               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 +866 @@
 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 +890 @@
      *                               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)

LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/top_bound_p.F

@@ -1 @@
-      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 +6 @@
 #include "paramet.h"
 #include "comconst.h"
-CC#include "comgeom2.h"
+#include "comvert.h"
+#include "comgeom2.h"
 
 
@@ -28 +29 @@
 c   -------------
 
-#include "comgeom.h"
 #include "comdissipn.h"
 
@@ -35 +35 @@
 
       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 +46 @@
       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 +86 @@
       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 +119 @@
         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 +149 @@
         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