source: trunk/libf/dyn3d/top_bound.F @ 108

Last change on this file since 108 was 108, checked in by slebonnois, 14 years ago

Sebastien Lebonnois: sponge layer et dissip horizontale.

File size: 4.2 KB
Line 
1      SUBROUTINE top_bound( vcov,ucov,teta,masse, du,dv,dh )
2      IMPLICIT NONE
3c
4#include "dimensions.h"
5#include "paramet.h"
6#include "comconst.h"
7#include "comvert.h"
8#include "comgeom2.h"
9
10
11c ..  DISSIPATION LINEAIRE A HAUT NIVEAU, RUN MESO,
12C     F. LOTT DEC. 2006
13c                                 (  10/12/06  )
14
15c=======================================================================
16c
17c   Auteur:  F. LOTT 
18c   -------
19c
20c   Objet:
21c   ------
22c
23c   Dissipation linéaire (ex top_bound de la physique)
24c
25c=======================================================================
26c-----------------------------------------------------------------------
27c   Declarations:
28c   -------------
29
30! #include "comgeom.h"
31#include "comdissipn.h"
32
33c   Arguments:
34c   ----------
35
36      REAL ucov(iip1,jjp1,llm),vcov(iip1,jjm,llm),teta(iip1,jjp1,llm)
37      REAL masse(iip1,jjp1,llm)
38      REAL dv(iip1,jjm,llm),du(iip1,jjp1,llm),dh(iip1,jjp1,llm)
39
40c   Local:
41c   ------
42
43      REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm),zm
44      REAL uzon(jjp1,llm),vzon(jjm,llm),tzon(jjp1,llm)
45     
46      INTEGER NDAMP
47      PARAMETER (NDAMP=4)
48      integer i
49      REAL,SAVE :: rdamp(llm)
50!     &   (/(0., i =1,llm-NDAMP),0.125E-5,.25E-5,.5E-5,1.E-5/)
51
52      LOGICAL,SAVE :: first=.true.
53
54      REAL zkm
55      INTEGER j,l
56
57
58C  CALCUL DES CHAMPS EN MOYENNE ZONALE:
59     
60      if (iflag_top_bound.eq.0) return
61
62      if (first) then
63         if (iflag_top_bound.eq.1) then
64! couche eponge dans les 4 dernieres couches du modele
65             rdamp(:)=0.
66             rdamp(llm)=tau_top_bound
67             rdamp(llm-1)=tau_top_bound/2.
68             rdamp(llm-2)=tau_top_bound/4.
69             rdamp(llm-3)=tau_top_bound/8.
70         else if (iflag_top_bound.eq.2) then
71! couche eponge dans toutes les couches de pression plus faible que
72! 100 fois la pression de la derniere couche
73             rdamp(:)=tau_top_bound
74     s       *max(presnivs(llm)/presnivs(:)-0.01,0.)
75         endif
76         first=.false.
77         print*,'TOP_BOUND mode',mode_top_bound
78         print*,'Coeffs pour la couche eponge a l equateur'
79         print*,'p (Pa)  z(km)  tau (s)   dt*rdamp'
80         do l=1,llm
81           if (rdamp(l).ne.0.) then
82            zkm        = phi(iip1/2,jjp1/2,l)/(1000*g)
83          print*,presnivs(l),zkm,
84     .          1./rdamp(l),
85     .          dt*rdamp(l)
86           endif
87         enddo
88      endif
89
90      CALL massbar(masse,massebx,masseby)
91
92c   mode = 0 : pas de sponge
93c   mode = 1 : u et v -> 0
94c   mode = 2 : u et v -> moyenne zonale
95c   mode = 3 : u, v et h -> moyenne zonale
96
97      if (mode_top_bound.ge.2) then
98       do l=1,llm
99        do j=1,jjm
100          vzon(j,l)=0.
101          zm=0.
102          do i=1,iim
103! Rm: on peut travailler directement avec la moyenne zonale de vcov
104! plutot qu'avec celle de v car le coefficient cv qui relie les deux
105! ne varie qu'en latitude
106            vzon(j,l)=vzon(j,l)+vcov(i,j,l)*masseby(i,j,l)
107            zm=zm+masseby(i,j,l)
108          enddo
109          vzon(j,l)=vzon(j,l)/zm
110        enddo
111       enddo
112
113       do l=1,llm
114        do j=2,jjm
115          uzon(j,l)=0.
116          zm=0.
117          do i=1,iim
118            uzon(j,l)=uzon(j,l)+massebx(i,j,l)*ucov(i,j,l)/cu(i,j)
119            zm=zm+massebx(i,j,l)
120          enddo
121          uzon(j,l)=uzon(j,l)/zm
122        enddo
123       enddo
124      else
125          vzon(:,:)=0.
126          uzon(:,:)=0.
127      endif
128
129      if (mode_top_bound.ge.3) then
130       do l=1,llm
131        do j=2,jjm
132          zm=0.
133          tzon(j,l)=0.
134          do i=1,iim
135            tzon(j,l)=tzon(j,l)+teta(i,j,l)*masse(i,j,l)
136            zm=zm+masse(i,j,l)
137          enddo
138          tzon(j,l)=tzon(j,l)/zm
139        enddo
140       enddo
141      endif
142
143C   AMORTISSEMENTS LINEAIRES:
144
145      if (mode_top_bound.ge.1) then
146       do l=1,llm
147        do i=1,iip1
148          do j=1,jjm
149            dv(i,j,l)= -rdamp(l)*(vcov(i,j,l)-vzon(j,l))
150          enddo
151        enddo
152       enddo
153
154       do l=1,llm
155        do i=1,iip1
156          do j=2,jjm
157            du(i,j,l)= -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
158          enddo
159        enddo
160       enddo
161      endif
162
163      if (mode_top_bound.ge.3) then
164       do l=1,llm
165        do i=1,iip1
166          do j=2,jjm
167            dh(i,j,l)= -rdamp(l)*(teta(i,j,l)-tzon(j,l))
168          enddo
169        enddo
170       enddo
171      endif
172     
173      RETURN
174      END
Note: See TracBrowser for help on using the repository browser.